# Load Necessary Packages
library(foreign)
library(stargazer)
library(Amelia)
library(Zelig)
library(dplyr)
library(MuMIn)
library(haven)
library(ggplot2)
library(grid)
library(gridExtra)
library(psych)
library(memisc)
library(ggpubr)

# Rescale Command
rescale <- function(x) {
  if(min(x,na.rm=T)>=0){
    x <- (x-min(x,na.rm=T))/(max(x,na.rm=T) - min(x,na.rm=T))
  } else{
    x <- x + abs(0 - min(x,na.rm=T))
    x <- (x-min(x,na.rm=T))/(max(x,na.rm=T) - min(x,na.rm=T))
  }
  print(x)
}

# Load Data
setwd("/users/josephphillips/Dropbox/Projects/Sorting and Affective Polarization/")
anes.cumulative <- read.dta("anes_timeseries_cdf.dta",convert.factors=T)
anes.cumulative <- subset(anes.cumulative,VCF0303=="1. Democrats (including leaners)" | VCF0303=="3. Republicans (including leaners)")

# Clean Feeling Thermometers
anes.cumulative$VCF0201[anes.cumulative$VCF0201>97] <- NA
anes.cumulative$VCF0202[anes.cumulative$VCF0202>97] <- NA
anes.cumulative$VCF0101[anes.cumulative$VCF0101==0] <- NA
anes.cumulative$VCF0316[anes.cumulative$VCF0316>5] <- NA
anes.cumulative$VCF0320[anes.cumulative$VCF0320>5] <- NA

# Create Affective Polarization Variable
anes.cumulative$polariz.old <- ifelse(as.numeric(anes.cumulative$VCF0303)==2,(anes.cumulative$VCF0201-anes.cumulative$VCF0202),ifelse(as.numeric(anes.cumulative$VCF0303)==4,(anes.cumulative$VCF0202-anes.cumulative$VCF0201),NA))
anes.cumulative$polariz.new <- ifelse(as.numeric(anes.cumulative$VCF0303)==2,(anes.cumulative$VCF0218-anes.cumulative$VCF0224),ifelse(as.numeric(anes.cumulative$VCF0303)==4,(anes.cumulative$VCF0224-anes.cumulative$VCF0218),NA))
anes.cumulative$polariz.likes <- ifelse(as.numeric(anes.cumulative$VCF0303)==2,(anes.cumulative$VCF0316-anes.cumulative$VCF0320),ifelse(as.numeric(anes.cumulative$VCF0303)==4,(anes.cumulative$VCF0320-anes.cumulative$VCF0316),NA))
anes.cumulative$polariz.combined <- ifelse(anes.cumulative$VCF0004<1978,anes.cumulative$polariz.old,anes.cumulative$polariz.new)
anes.cumulative$inparty.new <- ifelse(as.numeric(anes.cumulative$VCF0303)==2,(anes.cumulative$VCF0218),ifelse(as.numeric(anes.cumulative$VCF0303)==4,(anes.cumulative$VCF0224),NA))
anes.cumulative$outparty.new <- ifelse(as.numeric(anes.cumulative$VCF0303)==2,(anes.cumulative$VCF0224),ifelse(as.numeric(anes.cumulative$VCF0303)==4,(anes.cumulative$VCF0218),NA))
anes.cumulative$inparty.old <- ifelse(as.numeric(anes.cumulative$VCF0303)==2,(anes.cumulative$VCF0201),ifelse(as.numeric(anes.cumulative$VCF0303)==4,(anes.cumulative$VCF0202),NA))
anes.cumulative$outparty.old <- ifelse(as.numeric(anes.cumulative$VCF0303)==2,(anes.cumulative$VCF0202),ifelse(as.numeric(anes.cumulative$VCF0303)==4,(anes.cumulative$VCF0201),NA))
anes.cumulative$inparty.likes <- ifelse(as.numeric(anes.cumulative$VCF0303)==2,(anes.cumulative$VCF0316),ifelse(as.numeric(anes.cumulative$VCF0303)==4,(anes.cumulative$VCF0320),NA))
anes.cumulative$outparty.likes <- ifelse(as.numeric(anes.cumulative$VCF0303)==2,(anes.cumulative$VCF0320),ifelse(as.numeric(anes.cumulative$VCF0303)==4,(anes.cumulative$VCF0316),NA))
anes.cumulative$inparty.combined <- ifelse(anes.cumulative$VCF0004<1978,anes.cumulative$inparty.old,anes.cumulative$inparty.new)
anes.cumulative$outparty.combined <- ifelse(anes.cumulative$VCF0004<1978,anes.cumulative$outparty.old,anes.cumulative$outparty.new)

# Predictors of Affective Polarization
## PID
anes.cumulative$pid <- ifelse(anes.cumulative$VCF0301=="0. DK; NA; other; refused to answer; no Pre IW",NA,as.numeric(anes.cumulative$VCF0301) - 1)
## Strong Identifier
anes.cumulative$strong <- ifelse(anes.cumulative$pid==1 | anes.cumulative$pid==7,1,0)
## Age
anes.cumulative$age <- anes.cumulative$VCF0101-18
anes.cumulative$agesq <- anes.cumulative$age * anes.cumulative$age
## Gender
anes.cumulative$gender <- ifelse(anes.cumulative$VCF0104==0,NA,ifelse(anes.cumulative$VCF0104==8,0.5,ifelse(anes.cumulative$VCF0104==2,1,0)))
## Race
anes.cumulative$black <- ifelse(anes.cumulative$VCF0106=="0. Missing, pre-1966 data" | anes.cumulative$VCF0106=="9. Missing, DK/REF/NA",NA,ifelse(anes.cumulative$VCF0106=="2. Black non-Hispanic",1,0))
anes.cumulative$blacksort <- ifelse(anes.cumulative$pid<4 & anes.cumulative$black==1,1,ifelse(anes.cumulative$pid>4 & anes.cumulative$black==1,-1,0))
## Education
anes.cumulative$education <- ifelse(anes.cumulative$VCF0140=="8. DK" | anes.cumulative$VCF0140=="9. NA; RF; no Pre IW; short-form 'new' Cross Section",NA,as.numeric(anes.cumulative$VCF0140))
## South
anes.cumulative$south <- ifelse(anes.cumulative$VCF0112=="3. South (AL, AR, DE, D.C., FL, GA, KY, LA, MD, MS, NC",1,0)
## Ideological Sorting
anes.cumulative$ideo <- ifelse(anes.cumulative$VCF0803=="0. NA; no Post IW; form III,IV (1972); R not" | anes.cumulative$VCF0803=="9. DK; haven't thought much about it",NA,as.numeric(anes.cumulative$VCF0803)-1)
anes.cumulative$stronglib <- ifelse(anes.cumulative$VCF0803=="9. DK; haven't thought much about it" | anes.cumulative$VCF0803=="0. NA; no Post IW; form III,IV (1972); R not",NA,ifelse(anes.cumulative$VCF0803=="1. Extremely liberal",1,0))
anes.cumulative$strongcon <- ifelse(anes.cumulative$VCF0803=="9. DK; haven't thought much about it" | anes.cumulative$VCF0803=="0. NA; no Post IW; form III,IV (1972); R not",NA,ifelse(anes.cumulative$VCF0803=="7. Extremely conservative",1,0))
anes.cumulative$strongsort <- ifelse(anes.cumulative$pid<4 & anes.cumulative$stronglib==1,1,ifelse(anes.cumulative$pid<4 & anes.cumulative$strongcon==1,-1,ifelse(anes.cumulative$pid>4 & anes.cumulative$strongcon==1,1,ifelse(anes.cumulative$pid>4 & anes.cumulative$stronglib==1,-1,0))))
anes.cumulative$ideosort <- rescale(-abs(anes.cumulative$pid-anes.cumulative$ideo))
## Church Attendance
anes.cumulative$attendance <- as.numeric(anes.cumulative$VCF0130)
anes.cumulative$attendance[anes.cumulative$attendance==1 | anes.cumulative$attendance==6 | anes.cumulative$attendance==7] <- 1
anes.cumulative$attendance[anes.cumulative$attendance==5] <- -2
anes.cumulative$attendance[anes.cumulative$attendance==4] <- -3
anes.cumulative$attendance[anes.cumulative$attendance==3] <- -4
anes.cumulative$attendance[anes.cumulative$attendance==2] <- -5
anes.cumulative$attendance[anes.cumulative$attendance==-2] <- 2
anes.cumulative$attendance[anes.cumulative$attendance==-3] <- 3
anes.cumulative$attendance[anes.cumulative$attendance==-4] <- 4
anes.cumulative$attendance[anes.cumulative$attendance==-5] <- 5
anes.cumulative$attendance[anes.cumulative$attendance>7] <- NA
anes.cumulative$attendance <- (rescale(anes.cumulative$attendance)*2)-1
anes.cumulative$churchsort <- ifelse(anes.cumulative$pid>4,anes.cumulative$attendance,-anes.cumulative$attendance)
## Social Sorting
anes.cumulative$socialsort <- rescale(anes.cumulative$strong+anes.cumulative$blacksort+anes.cumulative$strongsort+anes.cumulative$churchsort)
## Interest (Media Proxy)
anes.cumulative$media <- ifelse(anes.cumulative$VCF0310=="0. NA; no Pre IW; version B (2008)" | anes.cumulative$VCF0310=="9. DK",NA,as.numeric(anes.cumulative$VCF0310)-2)
## Year
anes.cumulative$year <- anes.cumulative$VCF0004
anes.cumulative$cable <- ifelse(anes.cumulative$year>=1996,1,0)
anes.cumulative$year <- as.factor(anes.cumulative$VCF0004)

## APC Variables
anes.cumulative$agegroup <- ifelse(anes.cumulative$age<=6,20,ifelse(anes.cumulative$age>6 & anes.cumulative$age<=11,25,ifelse(anes.cumulative$age>11 & anes.cumulative$age<=16,30,ifelse(anes.cumulative$age>16 & anes.cumulative$age<=21,35,ifelse(anes.cumulative$age>21 & anes.cumulative$age<=26,40,ifelse(anes.cumulative$age>26 & anes.cumulative$age<=31,45,ifelse(anes.cumulative$age>31 & anes.cumulative$age<=36,50,ifelse(anes.cumulative$age>36 & anes.cumulative$age<=41,55,ifelse(anes.cumulative$age>41 & anes.cumulative$age<=46,60,ifelse(anes.cumulative$age>46 & anes.cumulative$age<=51,65,ifelse(anes.cumulative$age>51 & anes.cumulative$age<=56,70,ifelse(anes.cumulative$age>56 & anes.cumulative$age<=61,75,ifelse(anes.cumulative$age>61 & anes.cumulative$age<=66,80,ifelse(anes.cumulative$age>66 & anes.cumulative$age<=71,85,ifelse(anes.cumulative$age>71,90,NA)))))))))))))))
anes.cumulative$yeargroup <- ifelse(anes.cumulative$VCF0004>1975 & anes.cumulative$VCF0004<=1980,1980,ifelse(anes.cumulative$VCF0004>1980 & anes.cumulative$VCF0004<=1985,1985,ifelse(anes.cumulative$VCF0004>1985 & anes.cumulative$VCF0004<=1990,1990,ifelse(anes.cumulative$VCF0004>1990 & anes.cumulative$VCF0004<=1995,1995,ifelse(anes.cumulative$VCF0004>1995 & anes.cumulative$VCF0004<=2000,2000,ifelse(anes.cumulative$VCF0004>2000 & anes.cumulative$VCF0004<=2005,2005,ifelse(anes.cumulative$VCF0004>2005 & anes.cumulative$VCF0004<=2010,2010,ifelse(anes.cumulative$VCF0004>2010 & anes.cumulative$VCF0004<=2015,2015,ifelse(anes.cumulative$VCF0004>2015,2020,NA)))))))))
anes.cumulative$yeargroup.old <- ifelse(anes.cumulative$VCF0004>1960 & anes.cumulative$VCF0004<=1965,1965,ifelse(anes.cumulative$VCF0004>1965 & anes.cumulative$VCF0004<=1970,1970,ifelse(anes.cumulative$VCF0004>1970 & anes.cumulative$VCF0004<=1975,1975,ifelse(anes.cumulative$VCF0004>1975 & anes.cumulative$VCF0004<=1980,1980,ifelse(anes.cumulative$VCF0004>1980 & anes.cumulative$VCF0004<=1985,1985,ifelse(anes.cumulative$VCF0004>1985 & anes.cumulative$VCF0004<=1990,1990,ifelse(anes.cumulative$VCF0004>1990 & anes.cumulative$VCF0004<=1995,1995,ifelse(anes.cumulative$VCF0004>1995 & anes.cumulative$VCF0004<=2000,2000,ifelse(anes.cumulative$VCF0004>2000 & anes.cumulative$VCF0004<=2005,2005,ifelse(anes.cumulative$VCF0004>2005 & anes.cumulative$VCF0004<=2010,2010,ifelse(anes.cumulative$VCF0004>2010 & anes.cumulative$VCF0004<=2015,2015,ifelse(anes.cumulative$VCF0004>2015,2020,NA))))))))))))
anes.cumulative$yeargroup.likes <- ifelse(anes.cumulative$VCF0004>1950 & anes.cumulative$VCF0004<=1955,1955,ifelse(anes.cumulative$VCF0004>1955 & anes.cumulative$VCF0004<=1960,1960,ifelse(anes.cumulative$VCF0004>1960 & anes.cumulative$VCF0004<=1965,1965,ifelse(anes.cumulative$VCF0004>1965 & anes.cumulative$VCF0004<=1970,1970,ifelse(anes.cumulative$VCF0004>1970 & anes.cumulative$VCF0004<=1975,1975,ifelse(anes.cumulative$VCF0004>1975 & anes.cumulative$VCF0004<=1980,1980,ifelse(anes.cumulative$VCF0004>1980 & anes.cumulative$VCF0004<=1985,1985,ifelse(anes.cumulative$VCF0004>1985 & anes.cumulative$VCF0004<=1990,1990,ifelse(anes.cumulative$VCF0004>1990 & anes.cumulative$VCF0004<=1995,1995,ifelse(anes.cumulative$VCF0004>1995 & anes.cumulative$VCF0004<=2000,2000,ifelse(anes.cumulative$VCF0004>2000 & anes.cumulative$VCF0004<=2005,2005,ifelse(anes.cumulative$VCF0004>2005 & anes.cumulative$VCF0004<=2010,2010,ifelse(anes.cumulative$VCF0004>2010 & anes.cumulative$VCF0004<=2015,2015,ifelse(anes.cumulative$VCF0004>2015,2020,NA))))))))))))))
anes.cumulative$yeargroup.old.na <- ifelse(anes.cumulative$yeargroup.old>1985,NA,anes.cumulative$yeargroup.old)
anes.cumulative$yeargroup.likes.na <- ifelse(anes.cumulative$yeargroup.likes>2005,NA,anes.cumulative$yeargroup.likes)
anes.cumulative$cohort <- anes.cumulative$yeargroup-anes.cumulative$agegroup
anes.cumulative$cohort.old <- anes.cumulative$yeargroup.old.na-anes.cumulative$agegroup
anes.cumulative$cohort.combined <- anes.cumulative$yeargroup.old-anes.cumulative$agegroup
anes.cumulative$cohort.likes <- anes.cumulative$yeargroup.likes.na-anes.cumulative$agegroup
anes.cumulative$cohort.likes.full <- anes.cumulative$yeargroup.likes-anes.cumulative$agegroup
anes.cumulative$agegroup <- as.factor(anes.cumulative$agegroup)
anes.cumulative$yeargroup <- as.factor(anes.cumulative$yeargroup)
anes.cumulative$yeargroup.old <- as.factor(anes.cumulative$yeargroup.old)
anes.cumulative$yeargroup.likes <- as.factor(anes.cumulative$yeargroup.likes)
anes.cumulative$yeargroup.old.na <- as.factor(anes.cumulative$yeargroup.old.na)
anes.cumulative$yeargroup.likes.na <- as.factor(anes.cumulative$yeargroup.likes.na)
anes.cumulative$cohort <- as.factor(anes.cumulative$cohort)
anes.cumulative$cohort.old <- as.factor(anes.cumulative$cohort.old)
anes.cumulative$cohort.combined <- as.factor(anes.cumulative$cohort.combined)
anes.cumulative$cohort.likes <- as.factor(anes.cumulative$cohort.likes)

anes.cumulative.full <- anes.cumulative
anes.cumulative.dem <- subset(anes.cumulative.full,pid<=3)
anes.cumulative.rep <- subset(anes.cumulative.full,pid>=5)

# Table A3 (see models for age and period effects, cohort averages for cohort effects)
## Model 1
new.all.raw <- APCI::temp_model(data=anes.cumulative.full,outcome='polariz.new',age='agegroup',period='yeargroup',cohort='cohort',covariate=c(),family='gaussian')
new.all.raw.cohort <- APCI::cohortdeviation(new.all.raw$A,new.all.raw$P,new.all.raw$C,model=new.all.raw$model,covariate=c())
## Model 2
newpid.all.raw <- APCI::temp_model(data=anes.cumulative.full,outcome='polariz.new',age='agegroup',period='yeargroup',cohort='cohort',covariate=c('strong'),family='gaussian')
newpid.all.raw.cohort <- APCI::cohortdeviation(newpid.all.raw$A,newpid.all.raw$P,newpid.all.raw$C,model=newpid.all.raw$model,covariate=c('strong'))
## Model 3
newsort.all.raw <- APCI::temp_model(data=anes.cumulative.full,outcome='polariz.new',age='agegroup',period='yeargroup',cohort='cohort',covariate=c('ideosort','socialsort'),family='gaussian')
newsort.all.raw.cohort <- APCI::cohortdeviation(newsort.all.raw$A,newsort.all.raw$P,newsort.all.raw$C,model=newsort.all.raw$model,covariate=c('strong'))
## Model 4
new.all.controls <- APCI::temp_model(data=anes.cumulative.full,outcome='polariz.new',age='agegroup',period='yeargroup',cohort='cohort',covariate=c('gender','black','education','ideo','attendance','south'),family='gaussian')
new.all.controls.cohort <- APCI::cohortdeviation(new.all.controls$A,new.all.controls$P,new.all.controls$C,model=new.all.controls$model,covariate=c('gender','black','education','ideo','attendance','south'))
## Model 5
newpid.all.controls <- APCI::temp_model(data=anes.cumulative.full,outcome='polariz.new',age='agegroup',period='yeargroup',cohort='cohort',covariate=c('strong','gender','black','education','ideo','attendance','south'),family='gaussian')
newpid.all.controls.cohort <- APCI::cohortdeviation(newpid.all.controls$A,newpid.all.controls$P,newpid.all.controls$C,model=newpid.all.controls$model,covariate=c('strong','gender','black','education','ideo','attendance','south'))
## Model 6
newsort.all.controls <- APCI::temp_model(data=anes.cumulative.full,outcome='polariz.new',age='agegroup',period='yeargroup',cohort='cohort',covariate=c('ideosort','socialsort','gender','black','education','ideo','attendance','south'),family='gaussian')
newsort.all.controls.cohort <- APCI::cohortdeviation(newsort.all.controls$A,newsort.all.controls$P,newsort.all.controls$C,model=newsort.all.controls$model,covariate=c('ideosort','socialsort','gender','black','education','ideo','attendance','south'))

# Table A5 (see models for age and period effects, cohort averages for cohort effects)
## Model 1
new.dem.raw <- APCI::temp_model(data=anes.cumulative.dem,outcome='polariz.new',age='agegroup',period='yeargroup',cohort='cohort',covariate=c(),family='gaussian')
new.dem.raw.cohort <- APCI::cohortdeviation(new.dem.raw$A,new.dem.raw$P,new.dem.raw$C,model=new.dem.raw$model,covariate=c())
## Model 2
newpid.dem.raw <- APCI::temp_model(data=anes.cumulative.dem,outcome='polariz.new',age='agegroup',period='yeargroup',cohort='cohort',covariate=c('strong'),family='gaussian')
newpid.dem.raw.cohort <- APCI::cohortdeviation(newpid.dem.raw$A,newpid.dem.raw$P,newpid.dem.raw$C,model=newpid.dem.raw$model,covariate=c('strong'))
## Model 3
newsort.dem.raw <- APCI::temp_model(data=anes.cumulative.dem,outcome='polariz.new',age='agegroup',period='yeargroup',cohort='cohort',covariate=c('ideosort','socialsort'),family='gaussian')
newsort.dem.raw.cohort <- APCI::cohortdeviation(newsort.dem.raw$A,newsort.dem.raw$P,newsort.dem.raw$C,model=newsort.dem.raw$model,covariate=c('strong'))
## Model 4
new.dem.controls <- APCI::temp_model(data=anes.cumulative.dem,outcome='polariz.new',age='agegroup',period='yeargroup',cohort='cohort',covariate=c('gender','black','education','ideo','attendance','south'),family='gaussian')
new.dem.controls.cohort <- APCI::cohortdeviation(new.dem.controls$A,new.dem.controls$P,new.dem.controls$C,model=new.dem.controls$model,covariate=c('gender','black','education','ideo','attendance','south'))
## Model 5
newpid.dem.controls <- APCI::temp_model(data=anes.cumulative.dem,outcome='polariz.new',age='agegroup',period='yeargroup',cohort='cohort',covariate=c('strong','gender','black','education','ideo','attendance','south'),family='gaussian')
newpid.dem.controls.cohort <- APCI::cohortdeviation(newpid.dem.controls$A,newpid.dem.controls$P,newpid.dem.controls$C,model=newpid.dem.controls$model,covariate=c('strong','gender','black','education','ideo','attendance','south'))
## Model 6
newsort.dem.controls <- APCI::temp_model(data=anes.cumulative.dem,outcome='polariz.new',age='agegroup',period='yeargroup',cohort='cohort',covariate=c('ideosort','socialsort','gender','black','education','ideo','attendance','south'),family='gaussian')
newsort.dem.controls.cohort <- APCI::cohortdeviation(newsort.dem.controls$A,newsort.dem.controls$P,newsort.dem.controls$C,model=newsort.dem.controls$model,covariate=c('ideosort','socialsort','gender','black','education','ideo','attendance','south'))

# Table A7 (see models for age and period effects, cohort averages for cohort effects)
repnot90 <- subset(anes.cumulative.rep,agegroup!="90")
repnot90$agegroup <- as.factor(ifelse(repnot90$agegroup=="20","20",ifelse(repnot90$agegroup=="25","25",ifelse(repnot90$agegroup=="30","30",ifelse(repnot90$agegroup=="35","35",ifelse(repnot90$agegroup=="40","40",ifelse(repnot90$agegroup=="45","45",ifelse(repnot90$agegroup=="50","50",ifelse(repnot90$agegroup=="55","55",ifelse(repnot90$agegroup=="60","60",ifelse(repnot90$agegroup=="65","65",ifelse(repnot90$agegroup=="70","70",ifelse(repnot90$agegroup=="75","75",ifelse(repnot90$agegroup=="80","80","85"))))))))))))))
## Model 1
new.rep.raw <- APCI::temp_model(data=repnot90,outcome='polariz.new',age='agegroup',period='yeargroup',cohort='cohort',covariate=c(),family='gaussian')
new.rep.raw.cohort <- APCI::cohortdeviation(new.rep.raw$A,new.rep.raw$P,new.rep.raw$C,model=new.rep.raw$model,covariate=c())
## Model 2
newpid.rep.raw <- APCI::temp_model(data=repnot90,outcome='polariz.new',age='agegroup',period='yeargroup',cohort='cohort',covariate=c('strong'),family='gaussian')
newpid.rep.raw.cohort <- APCI::cohortdeviation(newpid.rep.raw$A,newpid.rep.raw$P,newpid.rep.raw$C,model=newpid.rep.raw$model,covariate=c('strong'))
## Model 3
newsort.rep.raw <- APCI::temp_model(data=repnot90,outcome='polariz.new',age='agegroup',period='yeargroup',cohort='cohort',covariate=c('ideosort','socialsort'),family='gaussian')
newsort.rep.raw.cohort <- APCI::cohortdeviation(newsort.rep.raw$A,newsort.rep.raw$P,newsort.rep.raw$C,model=newsort.rep.raw$model,covariate=c('strong'))
## Model 4
new.rep.controls <- APCI::temp_model(data=repnot90,outcome='polariz.new',age='agegroup',period='yeargroup',cohort='cohort',covariate=c('gender','black','education','ideo','attendance','south'),family='gaussian')
new.rep.controls.cohort <- APCI::cohortdeviation(new.rep.controls$A,new.rep.controls$P,new.rep.controls$C,model=new.rep.controls$model,covariate=c('gender','black','education','ideo','attendance','south'))
## Model 5
newpid.rep.controls <- APCI::temp_model(data=repnot90,outcome='polariz.new',age='agegroup',period='yeargroup',cohort='cohort',covariate=c('strong','gender','black','education','ideo','attendance','south'),family='gaussian')
newpid.rep.controls.cohort <- APCI::cohortdeviation(newpid.rep.controls$A,newpid.rep.controls$P,newpid.rep.controls$C,model=newpid.rep.controls$model,covariate=c('strong','gender','black','education','ideo','attendance','south'))
## Model 6
newsort.rep.controls <- APCI::temp_model(data=repnot90,outcome='polariz.new',age='agegroup',period='yeargroup',cohort='cohort',covariate=c('ideosort','socialsort','gender','black','education','ideo','attendance','south'),family='gaussian')
newsort.rep.controls.cohort <- APCI::cohortdeviation(newsort.rep.controls$A,newsort.rep.controls$P,newsort.rep.controls$C,model=newsort.rep.controls$model,covariate=c('ideosort','socialsort','gender','black','education','ideo','attendance','south'))

# Figure A1 (see Appendix Table A4 in Appendices Intrinsic Estimator.do for IE analysis coefficients)
appendixagenewie.dem <- data.frame(age=c(20,25,30,35,40,45,50,55,60,65,70,75,80,85,90),
                                   effect=c(-6.819,-5.951,-5.021,-3.522,-2.065,0.738,-0.825,1.801,1.398,3.134,3.874,0.464,0.993,2.587,9.214),
                                   se=c(1.285,1.134,1.012,0.918,0.890,0.878,0.66,0.895,0.969,1.054,1.235,1.444,1.727,2.358,4.544))
appendixagenewie.dem$lci <- appendixagenewie.dem$effect - (1.96*appendixagenewie.dem$se)
appendixagenewie.dem$uci <- appendixagenewie.dem$effect + (1.96*appendixagenewie.dem$se)
appendixagenewie.dem$color <- ifelse(appendixagenewie.dem$lci<0 & appendixagenewie.dem$uci>0,"Gray50","Black")
appendixagenewie.dem$group <- 1

appendixagenewapci.dem <- data.frame(age=c(25,30,35,40,45,50,55,60,65,70,75,80,85,90),
                                     effect=c(summary(new.dem.raw$model)$coefficients[2:15,1]),
                                     se=c(summary(new.dem.raw$model)$coefficients[2:15,2]))
appendixagenewapci.dem$lci <- appendixagenewapci.dem$effect - (1.96*appendixagenewapci.dem$se)
appendixagenewapci.dem$uci <- appendixagenewapci.dem$effect + (1.96*appendixagenewapci.dem$se)
appendixagenewapci.dem$color <- ifelse(appendixagenewapci.dem$lci<0 & appendixagenewapci.dem$uci>0,"Gray50","Black")
appendixagenewapci.dem$group <- 1

figurea1.a <- ggplot(appendixagenewapci.dem,aes(x=age,y=effect,color=color,group=group)) + geom_point() + geom_errorbar(aes(ymin=lci,ymax=uci),size=0.5,width=0) + theme_classic() + theme(plot.title = element_text(hjust = 0.5)) + ggtitle("APCI") + geom_hline(yintercept=0,linetype="dashed",color="red") + scale_color_identity() + xlab("Age") + ylab("Deviation from Grand Mean") + ylim(-25,52)
figurea1.b <- ggplot(appendixagenewie.dem,aes(x=age,y=effect,color=color,group=group)) + geom_point() + geom_errorbar(aes(ymin=lci,ymax=uci),size=0.5,width=0) + theme_classic() + theme(plot.title = element_text(hjust = 0.5)) + ggtitle("IE") + geom_hline(yintercept=0,linetype="dashed",color="red") + scale_color_identity() + xlab("Age") + ylab("Deviation from Grand Mean") + ylim(-25,52)
ggarrange(figurea1.b,figurea1.a,nrow=1)

# Figure A2
appendixperiodnewie.dem <- data.frame(period=c(1975,1980,1985,1990,1995,2000,2005,2010,2015),
                                      effect=c(-9.235,-3.570,-5.231,-5.933,-2.125,0.950,6.597,10.911,7.123),
                                      se=c(0.884,0.832,0.648,0.667,0.6,1.129,0.830,0.740,0.960))
appendixperiodnewie.dem$lci <- appendixperiodnewie.dem$effect - (1.96*appendixperiodnewie.dem$se)
appendixperiodnewie.dem$uci <- appendixperiodnewie.dem$effect + (1.96*appendixperiodnewie.dem$se)
appendixperiodnewie.dem$color <- ifelse(appendixperiodnewie.dem$lci<0 & appendixperiodnewie.dem$uci>0,"Gray50","Black")
appendixperiodnewie.dem$group <- 1

appendixperiodnewapci.dem <- data.frame(period=c(1980,1985,1990,1995,2000,2005,2010,2015),
                                        effect=c(summary(new.dem.raw$model)$coefficients[16:23,1]),
                                        se=c(summary(new.dem.raw$model)$coefficients[16:23,2]))
appendixperiodnewapci.dem$lci <- appendixperiodnewapci.dem$effect - (1.96*appendixperiodnewapci.dem$se)
appendixperiodnewapci.dem$uci <- appendixperiodnewapci.dem$effect + (1.96*appendixperiodnewapci.dem$se)
appendixperiodnewapci.dem$color <- ifelse(appendixperiodnewapci.dem$lci<0 & appendixperiodnewapci.dem$uci>0,"Gray50","Black")
appendixperiodnewapci.dem$group <- 1

figurea2.a <- ggplot(appendixperiodnewapci.dem,aes(x=period,y=effect,color=color,group=group)) + geom_point() + geom_errorbar(aes(ymin=lci,ymax=uci),size=0.5,width=0) + theme_classic() + theme(plot.title = element_text(hjust = 0.5)) + ggtitle("APCI") + geom_hline(yintercept=0,linetype="dashed",color="red") + scale_color_identity() + xlab("Period") + ylab("Deviation from Grand Mean") + ylim(-25,52)
figurea2.b <- ggplot(appendixperiodnewie.dem,aes(x=period,y=effect,color=color,group=group)) + geom_point() + geom_errorbar(aes(ymin=lci,ymax=uci),size=0.5,width=0) + theme_classic() + theme(plot.title = element_text(hjust = 0.5)) + ggtitle("IE") + geom_hline(yintercept=0,linetype="dashed",color="red") + scale_color_identity() + xlab("Period") + ylab("Deviation from Grand Mean") + ylim(-25,52)
ggarrange(figurea2.a,figurea2.b,nrow=1)

# Figure A3
appendixcohortnewie.dem <- data.frame(cohort=c(1885,1890,1895,1900,1905,1910,1915,1920,1925,1930,1935,1940,1945,1950,1955,1960,1965,1970,1975,1980,1985,1990,1995),
                                      effect=c(15.560,-7.101,2.258,5.562,0.607,-2.55,-0.373,-1.847,-1.646,-3.421,-5.616,-3.139,-2.918,-1.748,-0.853,0.056,0.949,-0.662,2.18,1.282,2.796,3.365,-0.843),
                                      se=c(12.176,8.974,4.642,3.460,2.608,2.282,2.06,1.856,1.685,1.586,1.425,1.234,1.044,0.887,0.776,0.764,0.845,0.998,1.113,1.256,1.411,1.617,2.66))
appendixcohortnewie.dem$lci <- appendixcohortnewie.dem$effect - (1.96*appendixcohortnewie.dem$se)
appendixcohortnewie.dem$uci <- appendixcohortnewie.dem$effect + (1.96*appendixcohortnewie.dem$se)
appendixcohortnewie.dem$color <- ifelse(appendixcohortnewie.dem$lci<0 & appendixcohortnewie.dem$uci>0,"Gray50","Black")
appendixcohortnewie.dem$group <- 1

appendixcohortnewapci.dem <- data.frame(cohort=c(1885,1890,1895,1900,1905,1910,1915,1920,1925,1930,1935,1940,1945,1950,1955,1960,1965,1970,1975,1980,1985,1990,1995),
                                        effect=c(as.numeric(new.dem.raw.cohort$cohort_average$cohort_average)),
                                        se=c(as.numeric(new.dem.raw.cohort$cohort_average$cohort_average_se)))
appendixcohortnewapci.dem$lci <- appendixcohortnewapci.dem$effect - (1.96*appendixcohortnewapci.dem$se)
appendixcohortnewapci.dem$uci <- appendixcohortnewapci.dem$effect + (1.96*appendixcohortnewapci.dem$se)
appendixcohortnewapci.dem$color <- ifelse(appendixcohortnewapci.dem$lci<0 & appendixcohortnewapci.dem$uci>0,"Gray50","Black")
appendixcohortnewapci.dem$group <- 1

figurea3.a <- ggplot(appendixcohortnewapci.dem,aes(x=cohort,y=effect,color=color,group=group)) + geom_point() + geom_errorbar(aes(ymin=lci,ymax=uci),size=0.5,width=0) + theme_classic() + theme(plot.title = element_text(hjust = 0.5)) + ggtitle("APCI") + geom_hline(yintercept=0,linetype="dashed",color="red") + scale_color_identity() + xlab("Birth Cohort") + ylab("Deviation from Grand Mean") + ylim(-25,52)
figurea3.b <- ggplot(appendixcohortnewie.dem,aes(x=cohort,y=effect,color=color,group=group)) + geom_point() + geom_errorbar(aes(ymin=lci,ymax=uci),size=0.5,width=0) + theme_classic() + theme(plot.title = element_text(hjust = 0.5)) + ggtitle("IE") + geom_hline(yintercept=0,linetype="dashed",color="red") + scale_color_identity() + xlab("Birth Cohort") + ylab("Deviation from Grand Mean") + ylim(-25,52)
ggarrange(figurea3.a,figurea3.b,nrow=1)

# Figure A4
appendixcohortnewapcislopes.dem <- data.frame(cohort=c(1885,1890,1895,1900,1905,1910,1915,1920,1925,1930,1935,1940,1945,1950,1955,1960,1965,1970,1975,1980,1985,1990,1995),
                                              effect=c(as.numeric(new.dem.raw.cohort$cohort_slope$cohort_slope)),
                                              se=c(as.numeric(new.dem.raw.cohort$cohort_slope$cohort_slope_se)))
appendixcohortnewapcislopes.dem$lci <- appendixcohortnewapcislopes.dem$effect - (1.96*appendixcohortnewapcislopes.dem$se)
appendixcohortnewapcislopes.dem$uci <- appendixcohortnewapcislopes.dem$effect + (1.96*appendixcohortnewapcislopes.dem$se)
appendixcohortnewapcislopes.dem$color <- ifelse(appendixcohortnewapcislopes.dem$lci<0 & appendixcohortnewapcislopes.dem$uci>0,"Gray50","Black")
appendixcohortnewapcislopes.dem$group <- 1

ggplot(appendixcohortnewie.dem,aes(x=cohort,y=effect,color=color,group=group)) + geom_point() + geom_errorbar(aes(ymin=lci,ymax=uci),size=0.5,width=0) + theme_classic() + theme(plot.title = element_text(hjust = 0.5)) + geom_hline(yintercept=0,linetype="dashed",color="red") + scale_color_identity() + xlab("Birth Cohort") + ylab("Deviation from Grand Mean") + ylim(-25,52)

# Figure A5
appendixagenewie.rep <- data.frame(age=c(20,25,30,35,40,45,50,55,60,65,70,75,80,85,90),
                                   effect=c(-5.894,-4.528,-5.323,-3.675,-1.584,-3.277,-2.859,-0.05,1.720,2.939,4.5,4.879,4.617,2.083,6.453),
                                   se=c(1.434,1.273,1.125,1.019,0.988,0.969,0.939,0.957,1.049,1.15,1.323,1.539,1.852,2.494,4.147))
appendixagenewie.rep$lci <- appendixagenewie.rep$effect - (1.96*appendixagenewie.rep$se)
appendixagenewie.rep$uci <- appendixagenewie.rep$effect + (1.96*appendixagenewie.rep$se)
appendixagenewie.rep$color <- ifelse(appendixagenewie.rep$lci<0 & appendixagenewie.rep$uci>0,"Gray50","Black")
appendixagenewie.rep$group <- 1

appendixagenewapci.rep <- data.frame(age=c(25,30,35,40,45,50,55,60,65,70,75,80,85),
                                     effect=c(summary(new.rep.raw$model)$coefficients[2:14,1]),
                                     se=c(summary(new.rep.raw$model)$coefficients[2:14,2]))
appendixagenewapci.rep$lci <- appendixagenewapci.rep$effect - (1.96*appendixagenewapci.rep$se)
appendixagenewapci.rep$uci <- appendixagenewapci.rep$effect + (1.96*appendixagenewapci.rep$se)
appendixagenewapci.rep$color <- ifelse(appendixagenewapci.rep$lci<0 & appendixagenewapci.rep$uci>0,"Gray50","Black")
appendixagenewapci.rep$group <- 1

figurea5.a <- ggplot(appendixagenewapci.rep,aes(x=age,y=effect,color=color,group=group)) + geom_point() + geom_errorbar(aes(ymin=lci,ymax=uci),size=0.5,width=0) + theme_classic() + theme(plot.title = element_text(hjust = 0.5)) + ggtitle("APCI") + geom_hline(yintercept=0,linetype="dashed",color="red") + scale_color_identity() + xlab("Age") + ylab("Deviation from Grand Mean") + ylim(-25,52) + xlim(20,90)
figurea5.b <- ggplot(appendixagenewie.rep,aes(x=age,y=effect,color=color,group=group)) + geom_point() + geom_errorbar(aes(ymin=lci,ymax=uci),size=0.5,width=0) + theme_classic() + theme(plot.title = element_text(hjust = 0.5)) + ggtitle("IE") + geom_hline(yintercept=0,linetype="dashed",color="red") + scale_color_identity() + xlab("Age") + ylab("Deviation from Grand Mean") + ylim(-25,52)
ggarrange(figurea5.a,figurea5.b,nrow=1)

# Figure A6
appendixperiodnewie.rep <- data.frame(period=c(1975,1980,1985,1990,1995,2000,2005,2010,2015),
                                      effect=c(-5.348,-2.417,-3.371,-1.383,-1.407,2.373,-2.515,6.983,7.085),
                                      se=c(1.011,0.921,0.714,0.713,0.666,1.191,1.087,0.828,0.996))
appendixperiodnewie.rep$lci <- appendixperiodnewie.rep$effect - (1.96*appendixperiodnewie.rep$se)
appendixperiodnewie.rep$uci <- appendixperiodnewie.rep$effect + (1.96*appendixperiodnewie.rep$se)
appendixperiodnewie.rep$color <- ifelse(appendixperiodnewie.rep$lci<0 & appendixperiodnewie.rep$uci>0,"Gray50","Black")
appendixperiodnewie.rep$group <- 1

appendixperiodnewapci.rep <- data.frame(period=c(1980,1985,1990,1995,2000,2005,2010,2015),
                                        effect=c(summary(new.rep.raw$model)$coefficients[15:22,1]),
                                        se=c(summary(new.rep.raw$model)$coefficients[15:22,2]))
appendixperiodnewapci.rep$lci <- appendixperiodnewapci.rep$effect - (1.96*appendixperiodnewapci.rep$se)
appendixperiodnewapci.rep$uci <- appendixperiodnewapci.rep$effect + (1.96*appendixperiodnewapci.rep$se)
appendixperiodnewapci.rep$color <- ifelse(appendixperiodnewapci.rep$lci<0 & appendixperiodnewapci.rep$uci>0,"Gray50","Black")
appendixperiodnewapci.rep$group <- 1

figurea6.a <- ggplot(appendixperiodnewapci.rep,aes(x=period,y=effect,color=color,group=group)) + geom_point() + geom_errorbar(aes(ymin=lci,ymax=uci),size=0.5,width=0) + theme_classic() + theme(plot.title = element_text(hjust = 0.5)) + ggtitle("APCI") + geom_hline(yintercept=0,linetype="dashed",color="red") + scale_color_identity() + xlab("Year") + ylab("Deviation from Grand Mean") + ylim(-25,52)
figurea6.b <- ggplot(appendixperiodnewie.rep,aes(x=period,y=effect,color=color,group=group)) + geom_point() + geom_errorbar(aes(ymin=lci,ymax=uci),size=0.5,width=0) + theme_classic() + theme(plot.title = element_text(hjust = 0.5)) + ggtitle("IE") + geom_hline(yintercept=0,linetype="dashed",color="red") + scale_color_identity() + xlab("Year") + ylab("Deviation from Grand Mean") + ylim(-25,52)
ggarrange(figurea6.a,figurea6.b,nrow=1)

# Figure A7
appendixcohortnewie.rep <- data.frame(cohort=c(1885,1890,1895,1900,1905,1910,1915,1920,1925,1930,1935,1940,1945,1950,1955,1960,1965,1970,1975,1980,1985,1990,1995),
                                      effect=c(23.463,-2.773,-3.481,-4.569,-1.886,-5.105,-3.646,-6.346,-4.244,1.612,0.256,-0.392,1.078,0.772,3.06,0.688,0.192,-0.312,-1.014,-0.151,-3.171,0.659,-1.652),
                                      se=c(14.503,7.532,4.373,3.680,2.939,2.521,2.309,2.101,1.836,1.705,1.523,1.324,1.145,0.996,0.857,0.821,0.917,1.105,1.341,1.614,1.659,1.922,3.006))
appendixcohortnewie.rep$lci <- appendixcohortnewie.rep$effect - (1.96*appendixcohortnewie.rep$se)
appendixcohortnewie.rep$uci <- appendixcohortnewie.rep$effect + (1.96*appendixcohortnewie.rep$se)
appendixcohortnewie.rep$color <- ifelse(appendixcohortnewie.rep$lci<0 & appendixcohortnewie.rep$uci>0,"Gray50","Black")
appendixcohortnewie.rep$group <- 1

appendixcohortnewapci.rep <- data.frame(cohort=c(1890,1895,1900,1905,1910,1915,1920,1925,1930,1935,1940,1945,1950,1955,1960,1965,1970,1975,1980,1985,1990,1995),
                                        effect=c(as.numeric(new.rep.raw.cohort$cohort_average$cohort_average)),
                                        se=c(as.numeric(new.rep.raw.cohort$cohort_average$cohort_average_se)))
appendixcohortnewapci.rep$lci <- appendixcohortnewapci.rep$effect - (1.96*appendixcohortnewapci.rep$se)
appendixcohortnewapci.rep$uci <- appendixcohortnewapci.rep$effect + (1.96*appendixcohortnewapci.rep$se)
appendixcohortnewapci.rep$color <- ifelse(appendixcohortnewapci.rep$lci<0 & appendixcohortnewapci.rep$uci>0,"Gray50","Black")
appendixcohortnewapci.rep$group <- 1

figurea7.a <- ggplot(appendixcohortnewapci.rep,aes(x=cohort,y=effect,color=color,group=group)) + geom_point() + geom_errorbar(aes(ymin=lci,ymax=uci),size=0.5,width=0) + theme_classic() + theme(plot.title = element_text(hjust = 0.5)) + ggtitle("APCI") + geom_hline(yintercept=0,linetype="dashed",color="red") + scale_color_identity() + xlab("Birth Cohort") + ylab("Deviation from Grand Mean") + ylim(-25,52)
figurea7.b <- ggplot(appendixcohortnewie.rep,aes(x=cohort,y=effect,color=color,group=group)) + geom_point() + geom_errorbar(aes(ymin=lci,ymax=uci),size=0.5,width=0) + theme_classic() + theme(plot.title = element_text(hjust = 0.5)) + ggtitle("IE") + geom_hline(yintercept=0,linetype="dashed",color="red") + scale_color_identity() + xlab("Birth Cohort") + ylab("Deviation from Grand Mean") + ylim(-25,52)
ggarrange(figurea7.a,figurea7.b,nrow=1)

# Figure A8
appendixcohortnewapcislopes.rep <- data.frame(cohort=c(1890,1895,1900,1905,1910,1915,1920,1925,1930,1935,1940,1945,1950,1955,1960,1965,1970,1975,1980,1985,1990,1995),
                                              effect=c(as.numeric(new.rep.raw.cohort$cohort_slope$cohort_slope)),
                                              se=c(as.numeric(new.rep.raw.cohort$cohort_slope$cohort_slope_se)))
appendixcohortnewapcislopes.rep$lci <- appendixcohortnewapcislopes.rep$effect - (1.96*appendixcohortnewapcislopes.rep$se)
appendixcohortnewapcislopes.rep$uci <- appendixcohortnewapcislopes.rep$effect + (1.96*appendixcohortnewapcislopes.rep$se)
appendixcohortnewapcislopes.rep$color <- ifelse(appendixcohortnewapcislopes.rep$lci<0 & appendixcohortnewapcislopes.rep$uci>0,"Gray50","Black")
appendixcohortnewapcislopes.rep$group <- 1

ggplot(appendixcohortnewie.rep,aes(x=cohort,y=effect,color=color,group=group)) + geom_point() + geom_errorbar(aes(ymin=lci,ymax=uci),size=0.5,width=0) + theme_classic() + theme(plot.title = element_text(hjust = 0.5)) + geom_hline(yintercept=0,linetype="dashed",color="red") + scale_color_identity() + xlab("Birth Cohort") + ylab("Deviation from Grand Mean") + ylim(-25,52)

# Table A9
anes.cumulative.full$yeargroup.sort <- anes.cumulative.full$yeargroup.old
anes.cumulative.full$yeargroup.sort[anes.cumulative.full$yeargroup.sort=="1965" | anes.cumulative.full$yeargroup.sort=="1970"] <- NA
anes.cumulative.full$yeargroup.sort <- droplevels(anes.cumulative.full$yeargroup.sort)
## Model 1
combined.all.raw <- APCI::temp_model(data=anes.cumulative.full,outcome='polariz.combined',age='agegroup',period='yeargroup.old',cohort='cohort.combined',covariate=c(),family='gaussian')
combined.all.raw.cohort <- APCI::cohortdeviation(combined.all.raw$A,combined.all.raw$P,combined.all.raw$C,model=combined.all.raw$model,covariate=c())
## Model 2
combinedpid.all.raw <- APCI::temp_model(data=anes.cumulative.full,outcome='polariz.combined',age='agegroup',period='yeargroup.old',cohort='cohort.combined',covariate=c('strong'),family='gaussian')
combinedpid.all.raw.cohort <- APCI::cohortdeviation(combinedpid.all.raw$A,combinedpid.all.raw$P,combinedpid.all.raw$C,model=combinedpid.all.raw$model,covariate=c('strong'))
## Model 3
combinedsort.all.raw <- APCI::temp_model(data=anes.cumulative.full,outcome='polariz.combined',age='agegroup',period='yeargroup.sort',cohort='cohort.combined',covariate=c('ideosort','socialsort'),family='gaussian')
combinedsort.all.raw.cohort <- APCI::cohortdeviation(combinedsort.all.raw$A,combinedsort.all.raw$P,combinedsort.all.raw$C,model=combinedsort.all.raw$model,covariate=c('strong'))
## Model 4
combined.all.controls <- APCI::temp_model(data=anes.cumulative.full,outcome='polariz.combined',age='agegroup',period='yeargroup.sort',cohort='cohort.combined',covariate=c('gender','black','education','ideo','attendance','south'),family='gaussian')
combined.all.controls.cohort <- APCI::cohortdeviation(combined.all.controls$A,combined.all.controls$P,combined.all.controls$C,model=combined.all.controls$model,covariate=c('gender','black','education','ideo','attendance','south'))
## Model 5
combinedpid.all.controls <- APCI::temp_model(data=anes.cumulative.full,outcome='polariz.combined',age='agegroup',period='yeargroup.sort',cohort='cohort.combined',covariate=c('strong','gender','black','education','ideo','attendance','south'),family='gaussian')
combinedpid.all.controls.cohort <- APCI::cohortdeviation(combinedpid.all.controls$A,combinedpid.all.controls$P,combinedpid.all.controls$C,model=combinedpid.all.controls$model,covariate=c('strong','gender','black','education','ideo','attendance','south'))
## Model 6
combinedsort.all.controls <- APCI::temp_model(data=anes.cumulative.full,outcome='polariz.combined',age='agegroup',period='yeargroup.sort',cohort='cohort.combined',covariate=c('ideosort','socialsort','gender','black','education','ideo','attendance','south'),family='gaussian')
combinedsort.all.controls.cohort <- APCI::cohortdeviation(combinedsort.all.controls$A,combinedsort.all.controls$P,combinedsort.all.controls$C,model=combinedsort.all.controls$model,covariate=c('ideosort','socialsort','gender','black','education','ideo','attendance','south'))

# Table A11
## Model 1
combined.dem.raw <- APCI::temp_model(data=anes.cumulative.dem,outcome='polariz.combined',age='agegroup',period='yeargroup.old',cohort='cohort.combined',covariate=c(),family='gaussian')
combined.dem.raw.cohort <- APCI::cohortdeviation(combined.dem.raw$A,combined.dem.raw$P,combined.dem.raw$C,model=combined.dem.raw$model,covariate=c())
## Model 2
combinedpid.dem.raw <- APCI::temp_model(data=anes.cumulative.dem,outcome='polariz.combined',age='agegroup',period='yeargroup.old',cohort='cohort.combined',covariate=c('strong'),family='gaussian')
combinedpid.dem.raw.cohort <- APCI::cohortdeviation(combinedpid.dem.raw$A,combinedpid.dem.raw$P,combinedpid.dem.raw$C,model=combinedpid.dem.raw$model,covariate=c('strong'))
## Model 3
combinedsort.dem.raw <- APCI::temp_model(data=demnot90,outcome='polariz.combined',age='agegroup',period='yeargroup.combined',cohort='cohort.combined',covariate=c('ideosort','socialsort'),family='gaussian')
combinedsort.dem.raw.cohort <- APCI::cohortdeviation(combinedsort.dem.raw$A,combinedsort.dem.raw$P,combinedsort.dem.raw$C,model=combinedsort.dem.raw$model,covariate=c('strong'))
## Model 4
combined.dem.controls <- APCI::temp_model(data=demnot90,outcome='polariz.combined',age='agegroup',period='yeargroup.combined',cohort='cohort.combined',covariate=c('gender','black','education','ideo','attendance','south'),family='gaussian')
combined.dem.controls.cohort <- APCI::cohortdeviation(combined.dem.controls$A,combined.dem.controls$P,combined.dem.controls$C,model=combined.dem.controls$model,covariate=c('gender','black','education','ideo','attendance','south'))
## Model 5
combinedpid.dem.controls <- APCI::temp_model(data=demnot90,outcome='polariz.combined',age='agegroup',period='yeargroup.combined',cohort='cohort.combined',covariate=c('strong','gender','black','education','ideo','attendance','south'),family='gaussian')
combinedpid.dem.controls.cohort <- APCI::cohortdeviation(combinedpid.dem.controls$A,combinedpid.dem.controls$P,combinedpid.dem.controls$C,model=combinedpid.dem.controls$model,covariate=c('strong','gender','black','education','ideo','attendance','south'))
## Model 6
combinedsort.dem.controls <- APCI::temp_model(data=demnot90,outcome='polariz.combined',age='agegroup',period='yeargroup.combined',cohort='cohort.combined',covariate=c('ideosort','socialsort','gender','black','education','ideo','attendance','south'),family='gaussian')
combinedsort.dem.controls.cohort <- APCI::cohortdeviation(combinedsort.dem.controls$A,combinedsort.dem.controls$P,combinedsort.dem.controls$C,model=combinedsort.dem.controls$model,covariate=c('ideosort','socialsort','gender','black','education','ideo','attendance','south'))

# Table A13
repnot90$yeargroup.combined1 <- repnot90$yeargroup.combined
repnot90$yeargroup.combined <- repnot90$yeargroup.old
repnot90$yeargroup.combined[repnot90$yeargroup.combined=="1965" | repnot90$yeargroup.combined=="1970"] <- NA
repnot90$yeargroup.combined <- droplevels(repnot90$yeargroup.combined)
## Model 1
combined.rep.raw <- APCI::temp_model(data=repnot90,outcome='polariz.combined',age='agegroup',period='yeargroup.old',cohort='cohort.combined',covariate=c(),family='gaussian')
combined.rep.raw.cohort <- APCI::cohortdeviation(combined.rep.raw$A,combined.rep.raw$P,combined.rep.raw$C,model=combined.rep.raw$model,covariate=c())
## Model 2
combinedpid.rep.raw <- APCI::temp_model(data=repnot90,outcome='polariz.combined',age='agegroup',period='yeargroup.old',cohort='cohort.combined',covariate=c('strong'),family='gaussian')
combinedpid.rep.raw.cohort <- APCI::cohortdeviation(combinedpid.rep.raw$A,combinedpid.rep.raw$P,combinedpid.rep.raw$C,model=combinedpid.rep.raw$model,covariate=c('strong'))
## Model 3
combinedsort.rep.raw <- APCI::temp_model(data=repnot90,outcome='polariz.combined',age='agegroup',period='yeargroup.combined',cohort='cohort.combined',covariate=c('ideosort','socialsort'),family='gaussian')
combinedsort.rep.raw.cohort <- APCI::cohortdeviation(combinedsort.rep.raw$A,combinedsort.rep.raw$P,combinedsort.rep.raw$C,model=combinedsort.rep.raw$model,covariate=c('strong'))
## Model 4
combined.rep.controls <- APCI::temp_model(data=repnot90,outcome='polariz.combined',age='agegroup',period='yeargroup.combined',cohort='cohort.combined',covariate=c('gender','black','education','ideo','attendance','south'),family='gaussian')
combined.rep.controls.cohort <- APCI::cohortdeviation(combined.rep.controls$A,combined.rep.controls$P,combined.rep.controls$C,model=combined.rep.controls$model,covariate=c('gender','black','education','ideo','attendance','south'))
## Model 5
combinedpid.rep.controls <- APCI::temp_model(data=repnot90,outcome='polariz.combined',age='agegroup',period='yeargroup.combined',cohort='cohort.combined',covariate=c('strong','gender','black','education','ideo','attendance','south'),family='gaussian')
combinedpid.rep.controls.cohort <- APCI::cohortdeviation(combinedpid.rep.controls$A,combinedpid.rep.controls$P,combinedpid.rep.controls$C,model=combinedpid.rep.controls$model,covariate=c('strong','gender','black','education','ideo','attendance','south'))
## Model 6
combinedsort.rep.controls <- APCI::temp_model(data=repnot90,outcome='polariz.combined',age='agegroup',period='yeargroup.combined',cohort='cohort.combined',covariate=c('ideosort','socialsort','gender','black','education','ideo','attendance','south'),family='gaussian')
combinedsort.rep.controls.cohort <- APCI::cohortdeviation(combinedsort.rep.controls$A,combinedsort.rep.controls$P,combinedsort.rep.controls$C,model=combinedsort.rep.controls$model,covariate=c('ideosort','socialsort','gender','black','education','ideo','attendance','south'))

# Figure A9
appendixagecombinedie.all <- data.frame(age=c(20,25,30,35,40,45,50,55,60,65,70,75,80,85,90),
                                        effect=c(-4.206,-3.528,-3.694,-2.404,-1.471,-0.746,-1.351,0.235,0.093,1.763,3.181,0.684,1.780,1.300,8.364),
                                        se=c(1.005,0.879,0.772,0.678,0.615,0.571,0.543,0.554,0.604,0.678,0.803,0.956,1.175,1.629,2.741))
appendixagecombinedie.all$lci <- appendixagecombinedie.all$effect - (1.96*appendixagecombinedie.all$se)
appendixagecombinedie.all$uci <- appendixagecombinedie.all$effect + (1.96*appendixagecombinedie.all$se)
appendixagecombinedie.all$color <- ifelse(appendixagecombinedie.all$lci<0 & appendixagecombinedie.all$uci>0,"Gray50","Black")
appendixagecombinedie.all$group <- 1

appendixagecombinedapci.all <- data.frame(age=c(25,30,35,40,45,50,55,60,65,70,75,80,85,90),
                                          effect=c(summary(combined.all.raw$model)$coefficients[2:15,1]),
                                          se=c(summary(combined.all.raw$model)$coefficients[2:15,2]))
appendixagecombinedapci.all$lci <- appendixagecombinedapci.all$effect - (1.96*appendixagecombinedapci.all$se)
appendixagecombinedapci.all$uci <- appendixagecombinedapci.all$effect + (1.96*appendixagecombinedapci.all$se)
appendixagecombinedapci.all$color <- ifelse(appendixagecombinedapci.all$lci<0 & appendixagecombinedapci.all$uci>0,"Gray50","Black")
appendixagecombinedapci.all$group <- 1

figurea9.a <- ggplot(appendixagecombinedapci.all,aes(x=age,y=effect,color=color,group=group)) + geom_point() + geom_errorbar(aes(ymin=lci,ymax=uci),size=0.5,width=0) + theme_classic() + theme(plot.title = element_text(hjust = 0.5)) + ggtitle("APCI") + geom_hline(yintercept=0,linetype="dashed",color="red") + scale_color_identity() + xlab("Age") + ylab("Deviation from Grand Mean") + ylim(-57,37)
figurea9.b <- ggplot(appendixagecombinedie.all,aes(x=age,y=effect,color=color,group=group)) + geom_point() + geom_errorbar(aes(ymin=lci,ymax=uci),size=0.5,width=0) + theme_classic() + theme(plot.title = element_text(hjust = 0.5)) + ggtitle("IE") + geom_hline(yintercept=0,linetype="dashed",color="red") + scale_color_identity() + xlab("Age") + ylab("Deviation from Grand Mean") + ylim(-57,37)
ggarrange(figurea9.a,figurea9.b,nrow=1)

# Figure A10
appendixperiodcombinedie.all <- data.frame(period=c(1960,1965,1970,1975,1980,1985,1990,1995,2000,2005,2010,2015),
                                           effect=c(0.273,-6.533,-13.241,-8.849,-1.332,-2.769,-2.092,0.348,3.977,7.092,12.817,10.309),
                                           se=c(0.970,0.750,0.689,0.546,0.547,0.414,0.471,0.481,0.88,0.764,0.721,0.877))
appendixperiodcombinedie.all$lci <- appendixperiodcombinedie.all$effect - (1.96*appendixperiodcombinedie.all$se)
appendixperiodcombinedie.all$uci <- appendixperiodcombinedie.all$effect + (1.96*appendixperiodcombinedie.all$se)
appendixperiodcombinedie.all$color <- ifelse(appendixperiodcombinedie.all$lci<0 & appendixperiodcombinedie.all$uci>0,"Gray50","Black")
appendixperiodcombinedie.all$group <- 1

appendixperiodcombinedapci.all <- data.frame(period=c(1965,1970,1975,1980,1985,1990,1995,2000,2005,2010,2015),
                                             effect=c(summary(combined.all.raw$model)$coefficients[16:26,1]),
                                             se=c(summary(combined.all.raw$model)$coefficients[16:26,2]))
appendixperiodcombinedapci.all$lci <- appendixperiodcombinedapci.all$effect - (1.96*appendixperiodcombinedapci.all$se)
appendixperiodcombinedapci.all$uci <- appendixperiodcombinedapci.all$effect + (1.96*appendixperiodcombinedapci.all$se)
appendixperiodcombinedapci.all$color <- ifelse(appendixperiodcombinedapci.all$lci<0 & appendixperiodcombinedapci.all$uci>0,"Gray50","Black")
appendixperiodcombinedapci.all$group <- 1

figurea10.a <- ggplot(appendixperiodcombinedapci.all,aes(x=period,y=effect,color=color,group=group)) + geom_point() + geom_errorbar(aes(ymin=lci,ymax=uci),size=0.5,width=0) + theme_classic() + theme(plot.title = element_text(hjust = 0.5)) + ggtitle("APCI") + geom_hline(yintercept=0,linetype="dashed",color="red") + scale_color_identity() + xlab("Year") + ylab("Deviation from Grand Mean") + ylim(-57,37)
figurea10.b <- ggplot(appendixperiodcombinedie.all,aes(x=period,y=effect,color=color,group=group)) + geom_point() + geom_errorbar(aes(ymin=lci,ymax=uci),size=0.5,width=0) + theme_classic() + theme(plot.title = element_text(hjust = 0.5)) + ggtitle("IE") + geom_hline(yintercept=0,linetype="dashed",color="red") + scale_color_identity() + xlab("Year") + ylab("Deviation from Grand Mean") + ylim(-57,37)
ggarrange(figurea10.a,figurea10.b,nrow=1)

# Figure A11
appendixcohortcombinedie.all <- data.frame(cohort=c(1870,1875,1880,1885,1890,1895,1900,1905,1910,1915,1920,1925,1930,1935,1940,1945,1950,1955,1960,1965,1970,1975,1980,1985,1990,1995),
                                           effect=c(2.505,-10.240,1.003,-0.534,-0.422,2.409,1.160,1.260,-0.053,0.655,-0.013,0.394,0.914,-0.548,0.205,0.393,0.903,1.474,0.627,-0.370,-0.852,0.756,0.521,-0.607,1.191,-2.732),
                                           se=c(12.690,8.648,4.947,2.987,2.456,2.091,1.894,1.687,1.520,1.379,1.244,1.11,1.007,0.890,0.761,0.645,0.561,0.516,0.531,0.621,0.761,0.896,1.041,1.149,1.320,2.034))
appendixcohortcombinedie.all$lci <- appendixcohortcombinedie.all$effect - (1.96*appendixcohortcombinedie.all$se)
appendixcohortcombinedie.all$uci <- appendixcohortcombinedie.all$effect + (1.96*appendixcohortcombinedie.all$se)
appendixcohortcombinedie.all$color <- ifelse(appendixcohortcombinedie.all$lci<0 & appendixcohortcombinedie.all$uci>0,"Gray50","Black")
appendixcohortcombinedie.all$group <- 1

appendixcohortcombinedapci.all <- data.frame(cohort=c(1870,1885,1880,1885,1890,1895,1900,1905,1910,1915,1920,1925,1930,1935,1940,1945,1950,1955,1960,1965,1970,1975,1980,1985,1990,1995),
                                             effect=c(as.numeric(combined.all.raw.cohort$cohort_average$cohort_average)),
                                             se=c(as.numeric(combined.all.raw.cohort$cohort_average$cohort_average_se)))
appendixcohortcombinedapci.all$lci <- appendixcohortcombinedapci.all$effect - (1.96*appendixcohortcombinedapci.all$se)
appendixcohortcombinedapci.all$uci <- appendixcohortcombinedapci.all$effect + (1.96*appendixcohortcombinedapci.all$se)
appendixcohortcombinedapci.all$color <- ifelse(appendixcohortcombinedapci.all$lci<0 & appendixcohortcombinedapci.all$uci>0,"Gray50","Black")
appendixcohortcombinedapci.all$group <- 1

figurea11.a <- ggplot(appendixcohortcombinedapci.all,aes(x=cohort,y=effect,color=color,group=group)) + geom_point() + geom_errorbar(aes(ymin=lci,ymax=uci),size=0.5,width=0) + theme_classic() + theme(plot.title = element_text(hjust = 0.5)) + ggtitle("APCI") + geom_hline(yintercept=0,linetype="dashed",color="red") + scale_color_identity() + xlab("Birth Cohort") + ylab("Deviation from Grand Mean") + ylim(-57,37)
figurea11.b <- ggplot(appendixcohortcombinedie.all,aes(x=cohort,y=effect,color=color,group=group)) + geom_point() + geom_errorbar(aes(ymin=lci,ymax=uci),size=0.5,width=0) + theme_classic() + theme(plot.title = element_text(hjust = 0.5)) + ggtitle("IE") + geom_hline(yintercept=0,linetype="dashed",color="red") + scale_color_identity() + xlab("Birth Cohort") + ylab("Deviation from Grand Mean") + ylim(-57,37)
ggarrange(figurea11.a,figurea11.b,nrow=1)

# Figure A12
appendixcohortcombinedapcislopes.all <- data.frame(cohort=c(1870,1875,1880,1885,1890,1895,1900,1905,1910,1915,1920,1925,1930,1935,1940,1945,1950,1955,1960,1965,1970,1975,1980,1985,1990,1995),
                                                   effect=c(as.numeric(combined.all.raw.cohort$cohort_slope$cohort_slope)),
                                                   se=c(as.numeric(combined.all.raw.cohort$cohort_slope$cohort_slope_se)))
appendixcohortcombinedapcislopes.all$lci <- appendixcohortcombinedapcislopes.all$effect - (1.96*appendixcohortcombinedapcislopes.all$se)
appendixcohortcombinedapcislopes.all$uci <- appendixcohortcombinedapcislopes.all$effect + (1.96*appendixcohortcombinedapcislopes.all$se)
appendixcohortcombinedapcislopes.all$color <- ifelse(appendixcohortcombinedapcislopes.all$lci<0 & appendixcohortcombinedapcislopes.all$uci>0,"Gray50","Black")
appendixcohortcombinedapcislopes.all$group <- 1

ggplot(appendixcohortcombinedapcislopes.all,aes(x=cohort,y=effect,color=color,group=group)) + geom_point() + geom_errorbar(aes(ymin=lci,ymax=uci),size=0.5,width=0) + theme_classic() + theme(plot.title = element_text(hjust = 0.5)) + geom_hline(yintercept=0,linetype="dashed",color="red") + scale_color_identity() + xlab("Birth Cohort") + ylab("Deviation from Grand Mean") + ylim(-57,37)

# Figure A13
appendixagecombinedie.dem <- data.frame(age=c(20,25,30,35,40,45,50,55,60,65,70,75,80,85,90),
                                        effect=c(-4.866,-4.511,-4.172,-2.734,-2.4,-0.37,-1.332,0.109,-0.116,1.387,3.695,-0.354,0.988,2.962,11.715),
                                        se=c(1.768,1.535,1.324,1.123,0.962,0.833,0.748,0.746,0.835,0.971,1.196,1.463,1.803,2.426,4.199))
appendixagecombinedie.dem$lci <- appendixagecombinedie.dem$effect - (1.96*appendixagecombinedie.dem$se)
appendixagecombinedie.dem$uci <- appendixagecombinedie.dem$effect + (1.96*appendixagecombinedie.dem$se)
appendixagecombinedie.dem$color <- ifelse(appendixagecombinedie.dem$lci<0 & appendixagecombinedie.dem$uci>0,"Gray50","Black")
appendixagecombinedie.dem$group <- 1

appendixagecombinedapci.dem <- data.frame(age=c(25,30,35,40,45,50,55,60,65,70,75,80,85,90),
                                          effect=c(summary(combined.dem.raw$model)$coefficients[2:15,1]),
                                          se=c(summary(combined.dem.raw$model)$coefficients[2:15,2]))
appendixagecombinedapci.dem$lci <- appendixagecombinedapci.dem$effect - (1.96*appendixagecombinedapci.dem$se)
appendixagecombinedapci.dem$uci <- appendixagecombinedapci.dem$effect + (1.96*appendixagecombinedapci.dem$se)
appendixagecombinedapci.dem$color <- ifelse(appendixagecombinedapci.dem$lci<0 & appendixagecombinedapci.dem$uci>0,"Gray50","Black")
appendixagecombinedapci.dem$group <- 1

figurea13.a <- ggplot(appendixagecombinedapci.dem,aes(x=age,y=effect,color=color,group=group)) + geom_point() + geom_errorbar(aes(ymin=lci,ymax=uci),size=0.5,width=0) + theme_classic() + theme(plot.title = element_text(hjust = 0.5)) + ggtitle("APCI") + geom_hline(yintercept=0,linetype="dashed",color="red") + scale_color_identity() + xlab("Age") + ylab("Deviation from Grand Mean") + ylim(-57,37)
figurea13.b <- ggplot(appendixagecombinedie.dem,aes(x=age,y=effect,color=color,group=group)) + geom_point() + geom_errorbar(aes(ymin=lci,ymax=uci),size=0.5,width=0) + theme_classic() + theme(plot.title = element_text(hjust = 0.5)) + ggtitle("IE") + geom_hline(yintercept=0,linetype="dashed",color="red") + scale_color_identity() + xlab("Age") + ylab("Deviation from Grand Mean") + ylim(-57,37)
ggarrange(figurea13.a,figurea13.b,nrow=1)

# Figure A14
appendixperiodcombinedie.dem <- data.frame(period=c(1960,1965,1970,1975,1980,1985,1990,1995,2000,2005,2010,2015),
                                           effect=c(0.951,-8.349,-12.728,-9.427,-1.185,-3.294,-3.873,0.157,3.482,9.459,14.149,10.659),
                                           se=c(1.554,1.24,1.083,0.830,0.769,0.561,0.647,0.69,1.273,1.138,1.217,1.503))
appendixperiodcombinedie.dem$lci <- appendixperiodcombinedie.dem$effect - (1.96*appendixperiodcombinedie.dem$se)
appendixperiodcombinedie.dem$uci <- appendixperiodcombinedie.dem$effect + (1.96*appendixperiodcombinedie.dem$se)
appendixperiodcombinedie.dem$color <- ifelse(appendixperiodcombinedie.dem$lci<0 & appendixperiodcombinedie.dem$uci>0,"Gray50","Black")
appendixperiodcombinedie.dem$group <- 1

appendixperiodcombinedapci.dem <- data.frame(period=c(1965,1970,1975,1980,1985,1990,1995,2000,2005,2010,2015),
                                             effect=c(summary(combined.dem.raw$model)$coefficients[16:26,1]),
                                             se=c(summary(combined.dem.raw$model)$coefficients[16:26,2]))
appendixperiodcombinedapci.dem$lci <- appendixperiodcombinedapci.dem$effect - (1.96*appendixperiodcombinedapci.dem$se)
appendixperiodcombinedapci.dem$uci <- appendixperiodcombinedapci.dem$effect + (1.96*appendixperiodcombinedapci.dem$se)
appendixperiodcombinedapci.dem$color <- ifelse(appendixperiodcombinedapci.dem$lci<0 & appendixperiodcombinedapci.dem$uci>0,"Gray50","Black")
appendixperiodcombinedapci.dem$group <- 1

figurea14.a <- ggplot(appendixperiodcombinedapci.dem,aes(x=period,y=effect,color=color,group=group)) + geom_point() + geom_errorbar(aes(ymin=lci,ymax=uci),size=0.5,width=0) + theme_classic() + theme(plot.title = element_text(hjust = 0.5)) + ggtitle("APCI") + geom_hline(yintercept=0,linetype="dashed",color="red") + scale_color_identity() + xlab("Year") + ylab("Deviation from Grand Mean") + ylim(-57,37)
figurea14.b <- ggplot(appendixperiodcombinedie.dem,aes(x=period,y=effect,color=color,group=group)) + geom_point() + geom_errorbar(aes(ymin=lci,ymax=uci),size=0.5,width=0) + theme_classic() + theme(plot.title = element_text(hjust = 0.5)) + ggtitle("IE") + geom_hline(yintercept=0,linetype="dashed",color="red") + scale_color_identity() + xlab("Year") + ylab("Deviation from Grand Mean") + ylim(-57,37)
ggarrange(figurea14.a,figurea14.b,nrow=1)

# Figure A15
appendixcohortcombinedie.dem <- data.frame(cohort=c(1870,1875,1880,1885,1890,1895,1900,1905,1910,1915,1920,1925,1930,1935,1940,1945,1950,1955,1960,1965,1970,1975,1980,1985,1990,1995),
                                           effect=c(0.310,0.699,0.048,-3.517,0.062,-0.390,2.287,-0.487,-0.399,0.842,0.015,0.496,-0.596,-2.279,-0.517,-0.631,0.253,0.123,0.815,-0.396,-0.359,1.939,0.672,1.819,1.888,-2.696),
                                           se=c(25.515,13.365,7.340,4.771,4.077,3.534,3.199,2.87,2.594,2.346,2.095,1.862,1.651,1.437,1.199,0.978,0.797,0.704,0.723,0.864,1.085,1.281,1.501,1.736,2.006,2.954))
appendixcohortcombinedie.dem$lci <- appendixcohortcombinedie.dem$effect - (1.96*appendixcohortcombinedie.dem$se)
appendixcohortcombinedie.dem$uci <- appendixcohortcombinedie.dem$effect + (1.96*appendixcohortcombinedie.dem$se)
appendixcohortcombinedie.dem$color <- ifelse(appendixcohortcombinedie.dem$lci<0 & appendixcohortcombinedie.dem$uci>0,"Gray50","Black")
appendixcohortcombinedie.dem$group <- 1

appendixcohortcombinedapci.dem <- data.frame(cohort=c(1870,1885,1880,1885,1890,1895,1900,1905,1910,1915,1920,1925,1930,1935,1940,1945,1950,1955,1960,1965,1970,1975,1980,1985,1990,1995),
                                             effect=c(as.numeric(combined.dem.raw.cohort$cohort_average$cohort_average)),
                                             se=c(as.numeric(combined.dem.raw.cohort$cohort_average$cohort_average_se)))
appendixcohortcombinedapci.dem$lci <- appendixcohortcombinedapci.dem$effect - (1.96*appendixcohortcombinedapci.dem$se)
appendixcohortcombinedapci.dem$uci <- appendixcohortcombinedapci.dem$effect + (1.96*appendixcohortcombinedapci.dem$se)
appendixcohortcombinedapci.dem$color <- ifelse(appendixcohortcombinedapci.dem$lci<0 & appendixcohortcombinedapci.dem$uci>0,"Gray50","Black")
appendixcohortcombinedapci.dem$group <- 1

figurea15.a <- ggplot(appendixcohortcombinedapci.dem,aes(x=cohort,y=effect,color=color,group=group)) + geom_point() + geom_errorbar(aes(ymin=lci,ymax=uci),size=0.5,width=0) + theme_classic() + theme(plot.title = element_text(hjust = 0.5)) + ggtitle("APCI") + geom_hline(yintercept=0,linetype="dashed",color="red") + scale_color_identity() + xlab("Birth Cohort") + ylab("Deviation from Grand Mean") + ylim(-57,37)
figurea15.b <- ggplot(appendixcohortcombinedie.dem,aes(x=cohort,y=effect,color=color,group=group)) + geom_point() + geom_errorbar(aes(ymin=lci,ymax=uci),size=0.5,width=0) + theme_classic() + theme(plot.title = element_text(hjust = 0.5)) + ggtitle("IE") + geom_hline(yintercept=0,linetype="dashed",color="red") + scale_color_identity() + xlab("Birth Cohort") + ylab("Deviation from Grand Mean") + ylim(-57,37)
ggarrange(figurea15.a,figurea15.b,nrow=1)

# Figure A16
appendixcohortcombinedapcislopes.dem <- data.frame(cohort=c(1870,1875,1880,1885,1890,1895,1900,1905,1910,1915,1920,1925,1930,1935,1940,1945,1950,1955,1960,1965,1970,1975,1980,1985,1990,1995),
                                                   effect=c(as.numeric(combined.dem.raw.cohort$cohort_slope$cohort_slope)),
                                                   se=c(as.numeric(combined.dem.raw.cohort$cohort_slope$cohort_slope_se)))
appendixcohortcombinedapcislopes.dem$lci <- appendixcohortcombinedapcislopes.dem$effect - (1.96*appendixcohortcombinedapcislopes.dem$se)
appendixcohortcombinedapcislopes.dem$uci <- appendixcohortcombinedapcislopes.dem$effect + (1.96*appendixcohortcombinedapcislopes.dem$se)
appendixcohortcombinedapcislopes.dem$color <- ifelse(appendixcohortcombinedapcislopes.dem$lci<0 & appendixcohortcombinedapcislopes.dem$uci>0,"Gray50","Black")
appendixcohortcombinedapcislopes.dem$group <- 1

ggplot(appendixcohortcombinedapcislopes.dem,aes(x=cohort,y=effect,color=color,group=group)) + geom_point() + geom_errorbar(aes(ymin=lci,ymax=uci),size=0.5,width=0) + theme_classic() + theme(plot.title = element_text(hjust = 0.5)) + geom_hline(yintercept=0,linetype="dashed",color="red") + scale_color_identity() + xlab("Birth Cohort") + ylab("Deviation from Grand Mean") + ylim(-57,37)

# Figure A17
appendixagecombinedie.rep <- data.frame(age=c(20,25,30,35,40,45,50,55,60,65,70,75,80,85,90),
                                        effect=c(-4.364,-3.204,-3.745,-2.671,-0.813,-1.932,-1.966,0.022,0.22,2.225,2.993,2.557,3.126,-0.654,8.204),
                                        se=c(1.284,1.126,1.001,0.903,0.849,0.821,0.796,0.820,0.877,0.974,1.125,1.305,1.601,2.265,3.59))
appendixagecombinedie.rep$lci <- appendixagecombinedie.rep$effect - (1.96*appendixagecombinedie.rep$se)
appendixagecombinedie.rep$uci <- appendixagecombinedie.rep$effect + (1.96*appendixagecombinedie.rep$se)
appendixagecombinedie.rep$color <- ifelse(appendixagecombinedie.rep$lci<0 & appendixagecombinedie.rep$uci>0,"Gray50","Black")
appendixagecombinedie.rep$group <- 1

appendixagecombinedapci.rep <- data.frame(age=c(25,30,35,40,45,50,55,60,65,70,75,80,85),
                                          effect=c(summary(combined.rep.raw$model)$coefficients[2:14,1]),
                                          se=c(summary(combined.rep.raw$model)$coefficients[2:14,2]))
appendixagecombinedapci.rep$lci <- appendixagecombinedapci.rep$effect - (1.96*appendixagecombinedapci.rep$se)
appendixagecombinedapci.rep$uci <- appendixagecombinedapci.rep$effect + (1.96*appendixagecombinedapci.rep$se)
appendixagecombinedapci.rep$color <- ifelse(appendixagecombinedapci.rep$lci<0 & appendixagecombinedapci.rep$uci>0,"Gray50","Black")
appendixagecombinedapci.rep$group <- 1

figurea17.a <- ggplot(appendixagecombinedapci.rep,aes(x=age,y=effect,color=color,group=group)) + geom_point() + geom_errorbar(aes(ymin=lci,ymax=uci),size=0.5,width=0) + theme_classic() + theme(plot.title = element_text(hjust = 0.5)) + ggtitle("APCI") + geom_hline(yintercept=0,linetype="dashed",color="red") + scale_color_identity() + xlab("Age") + ylab("Deviation from Grand Mean") + ylim(-57,37)
figurea17.b <- ggplot(appendixagecombinedie.rep,aes(x=age,y=effect,color=color,group=group)) + geom_point() + geom_errorbar(aes(ymin=lci,ymax=uci),size=0.5,width=0) + theme_classic() + theme(plot.title = element_text(hjust = 0.5)) + ggtitle("IE") + geom_hline(yintercept=0,linetype="dashed",color="red") + scale_color_identity() + xlab("Age") + ylab("Deviation from Grand Mean") + ylim(-57,37)
ggarrange(figurea17.a,figurea17.b,nrow=1)

# Figure A18
appendixperiodcombinedie.rep <- data.frame(period=c(1960,1965,1970,1975,1980,1985,1990,1995,2000,2005,2010,2015),
                                           effect=c(-1.395,-3.175,-13.385,-7.608,-1.014,-1.619,0.625,0.937,5.078,0.475,10.363,10.719),
                                           se=c(1.386,1.016,0.953,0.78,0.798,0.611,0.679,0.69,1.234,1.158,0.945,1.112))
appendixperiodcombinedie.rep$lci <- appendixperiodcombinedie.rep$effect - (1.96*appendixperiodcombinedie.rep$se)
appendixperiodcombinedie.rep$uci <- appendixperiodcombinedie.rep$effect + (1.96*appendixperiodcombinedie.rep$se)
appendixperiodcombinedie.rep$color <- ifelse(appendixperiodcombinedie.rep$lci<0 & appendixperiodcombinedie.rep$uci>0,"Gray50","Black")
appendixperiodcombinedie.rep$group <- 1

appendixperiodcombinedapci.rep <- data.frame(period=c(1965,1970,1975,1980,1985,1990,1995,2000,2005,2010,2015),
                                             effect=c(summary(combined.rep.raw$model)$coefficients[15:25,1]),
                                             se=c(summary(combined.rep.raw$model)$coefficients[15:25,2]))
appendixperiodcombinedapci.rep$lci <- appendixperiodcombinedapci.rep$effect - (1.96*appendixperiodcombinedapci.rep$se)
appendixperiodcombinedapci.rep$uci <- appendixperiodcombinedapci.rep$effect + (1.96*appendixperiodcombinedapci.rep$se)
appendixperiodcombinedapci.rep$color <- ifelse(appendixperiodcombinedapci.rep$lci<0 & appendixperiodcombinedapci.rep$uci>0,"Gray50","Black")
appendixperiodcombinedapci.rep$group <- 1

figurea18.a <- ggplot(appendixperiodcombinedapci.rep,aes(x=period,y=effect,color=color,group=group)) + geom_point() + geom_errorbar(aes(ymin=lci,ymax=uci),size=0.5,width=0) + theme_classic() + theme(plot.title = element_text(hjust = 0.5)) + ggtitle("APCI") + geom_hline(yintercept=0,linetype="dashed",color="red") + scale_color_identity() + xlab("Year") + ylab("Deviation from Grand Mean") + ylim(-57,37)
figurea18.b <- ggplot(appendixperiodcombinedie.rep,aes(x=period,y=effect,color=color,group=group)) + geom_point() + geom_errorbar(aes(ymin=lci,ymax=uci),size=0.5,width=0) + theme_classic() + theme(plot.title = element_text(hjust = 0.5)) + ggtitle("IE") + geom_hline(yintercept=0,linetype="dashed",color="red") + scale_color_identity() + xlab("Year") + ylab("Deviation from Grand Mean") + ylim(-57,37)
ggarrange(figurea18.a,figurea18.b,nrow=1)

# Figure A19
appendixcohortcombinedie.rep <- data.frame(cohort=c(1870,1875,1880,1885,1890,1895,1900,1905,1910,1915,1920,1925,1930,1935,1940,1945,1950,1955,1960,1965,1970,1975,1980,1985,1990,1995),
                                           effect=c(6.668,-19.587,2.403,2.484,-1.213,5.286,-1.216,2.613,-0.315,-0.693,-1.227,-0.249,3.093,1.934,1.398,1.992,1.774,3.717,1.076,0.260,-0.518,-1.323,-0.931,-4.156,-0.363,-2.905),
                                           se=c(14.215,11.143,6.79,3.956,3.159,2.69,2.469,2.209,1.979,1.8,1.651,1.459,1.352,1.191,1.038,0.908,0.826,0.756,0.771,0.896,1.094,1.335,1.601,1.651,1.908,2.933))
appendixcohortcombinedie.rep$lci <- appendixcohortcombinedie.rep$effect - (1.96*appendixcohortcombinedie.rep$se)
appendixcohortcombinedie.rep$uci <- appendixcohortcombinedie.rep$effect + (1.96*appendixcohortcombinedie.rep$se)
appendixcohortcombinedie.rep$color <- ifelse(appendixcohortcombinedie.rep$lci<0 & appendixcohortcombinedie.rep$uci>0,"Gray50","Black")
appendixcohortcombinedie.rep$group <- 1

appendixcohortcombinedapci.rep <- data.frame(cohort=c(1875,1880,1885,1890,1895,1900,1905,1910,1915,1920,1925,1930,1935,1940,1945,1950,1955,1960,1965,1970,1975,1980,1985,1990,1995),
                                             effect=c(as.numeric(combined.rep.raw.cohort$cohort_average$cohort_average)),
                                             se=c(as.numeric(combined.rep.raw.cohort$cohort_average$cohort_average_se)))
appendixcohortcombinedapci.rep$lci <- appendixcohortcombinedapci.rep$effect - (1.96*appendixcohortcombinedapci.rep$se)
appendixcohortcombinedapci.rep$uci <- appendixcohortcombinedapci.rep$effect + (1.96*appendixcohortcombinedapci.rep$se)
appendixcohortcombinedapci.rep$color <- ifelse(appendixcohortcombinedapci.rep$lci<0 & appendixcohortcombinedapci.rep$uci>0,"Gray50","Black")
appendixcohortcombinedapci.rep$group <- 1

figurea19.a <- ggplot(appendixcohortcombinedapci.rep,aes(x=cohort,y=effect,color=color,group=group)) + geom_point() + geom_errorbar(aes(ymin=lci,ymax=uci),size=0.5,width=0) + theme_classic() + theme(plot.title = element_text(hjust = 0.5)) + ggtitle("APCI") + geom_hline(yintercept=0,linetype="dashed",color="red") + scale_color_identity() + xlab("Birth Cohort") + ylab("Deviation from Grand Mean") + ylim(-57,37)
figurea19.b <- ggplot(appendixcohortcombinedie.rep,aes(x=cohort,y=effect,color=color,group=group)) + geom_point() + geom_errorbar(aes(ymin=lci,ymax=uci),size=0.5,width=0) + theme_classic() + theme(plot.title = element_text(hjust = 0.5)) + ggtitle("IE") + geom_hline(yintercept=0,linetype="dashed",color="red") + scale_color_identity() + xlab("Birth Cohort") + ylab("Deviation from Grand Mean") + ylim(-57,37)
ggarrange(figurea19.a,figurea19.b,nrow=1)

# Figure A20
appendixcohortcombinedapcislopes.rep <- data.frame(cohort=c(1875,1880,1885,1890,1895,1900,1905,1910,1915,1920,1925,1930,1935,1940,1945,1950,1955,1960,1965,1970,1975,1980,1985,1990,1995),
                                                   effect=c(as.numeric(combined.rep.raw.cohort$cohort_slope$cohort_slope)),
                                                   se=c(as.numeric(combined.rep.raw.cohort$cohort_slope$cohort_slope_se)))
appendixcohortcombinedapcislopes.rep$lci <- appendixcohortcombinedapcislopes.rep$effect - (1.96*appendixcohortcombinedapcislopes.rep$se)
appendixcohortcombinedapcislopes.rep$uci <- appendixcohortcombinedapcislopes.rep$effect + (1.96*appendixcohortcombinedapcislopes.rep$se)
appendixcohortcombinedapcislopes.rep$color <- ifelse(appendixcohortcombinedapcislopes.rep$lci<0 & appendixcohortcombinedapcislopes.rep$uci>0,"Gray50","Black")
appendixcohortcombinedapcislopes.rep$group <- 1

ggplot(appendixcohortcombinedapcislopes.rep,aes(x=cohort,y=effect,color=color,group=group)) + geom_point() + geom_errorbar(aes(ymin=lci,ymax=uci),size=0.5,width=0) + theme_classic() + theme(plot.title = element_text(hjust = 0.5)) + geom_hline(yintercept=0,linetype="dashed",color="red") + scale_color_identity() + xlab("Cohort") + ylab("Deviation from Grand Mean") + ylim(-57,37)

# Table A15
## Model 1
inpartynew.all.raw <- APCI::temp_model(data=anes.cumulative.full,outcome='inparty.new',age='agegroup',period='yeargroup',cohort='cohort',covariate=c(),family='gaussian')
inpartynew.all.raw.cohort <- APCI::cohortdeviation(inpartynew.all.raw$A,inpartynew.all.raw$P,inpartynew.all.raw$C,model=inpartynew.all.raw$model,covariate=c())
## Model 2
inpartynewpid.all.raw <- APCI::temp_model(data=anes.cumulative.full,outcome='inparty.new',age='agegroup',period='yeargroup',cohort='cohort',covariate=c('strong'),family='gaussian')
inpartynewpid.all.raw.cohort <- APCI::cohortdeviation(inpartynewpid.all.raw$A,inpartynewpid.all.raw$P,inpartynewpid.all.raw$C,model=inpartynewpid.all.raw$model,covariate=c('strong'))
## Model 3
inpartynewsort.all.raw <- APCI::temp_model(data=anes.cumulative.full,outcome='inparty.new',age='agegroup',period='yeargroup',cohort='cohort',covariate=c('ideosort','socialsort'),family='gaussian')
inpartynewsort.all.raw.cohort <- APCI::cohortdeviation(inpartynewsort.all.raw$A,inpartynewsort.all.raw$P,inpartynewsort.all.raw$C,model=inpartynewsort.all.raw$model,covariate=c('strong'))
## Model 4
inpartynew.all.controls <- APCI::temp_model(data=anes.cumulative.full,outcome='inparty.new',age='agegroup',period='yeargroup',cohort='cohort',covariate=c('gender','black','education','ideo','attendance','south'),family='gaussian')
inpartynew.all.controls.cohort <- APCI::cohortdeviation(inpartynew.all.controls$A,inpartynew.all.controls$P,inpartynew.all.controls$C,model=inpartynew.all.controls$model,covariate=c('gender','black','education','ideo','attendance','south'))
## Model 5
inpartynewpid.all.controls <- APCI::temp_model(data=anes.cumulative.full,outcome='inparty.new',age='agegroup',period='yeargroup',cohort='cohort',covariate=c('strong','gender','black','education','ideo','attendance','south'),family='gaussian')
inpartynewpid.all.controls.cohort <- APCI::cohortdeviation(inpartynewpid.all.controls$A,inpartynewpid.all.controls$P,inpartynewpid.all.controls$C,model=inpartynewpid.all.controls$model,covariate=c('strong','gender','black','education','ideo','attendance','south'))
## Model 6
inpartynewsort.all.controls <- APCI::temp_model(data=anes.cumulative.full,outcome='inparty.new',age='agegroup',period='yeargroup',cohort='cohort',covariate=c('ideosort','socialsort','gender','black','education','ideo','attendance','south'),family='gaussian')
inpartynewsort.all.controls.cohort <- APCI::cohortdeviation(inpartynewsort.all.controls$A,inpartynewsort.all.controls$P,inpartynewsort.all.controls$C,model=inpartynewsort.all.controls$model,covariate=c('ideosort','socialsort','gender','black','education','ideo','attendance','south'))

# Table A17
## Model 1
inpartynew.dem.raw <- APCI::temp_model(data=anes.cumulative.dem,outcome='inparty.new',age='agegroup',period='yeargroup',cohort='cohort',covariate=c(),family='gaussian')
inpartynew.dem.raw.cohort <- APCI::cohortdeviation(inpartynew.dem.raw$A,inpartynew.dem.raw$P,inpartynew.dem.raw$C,model=inpartynew.dem.raw$model,covariate=c())
## Model 2
inpartynewpid.dem.raw <- APCI::temp_model(data=anes.cumulative.dem,outcome='inparty.new',age='agegroup',period='yeargroup',cohort='cohort',covariate=c('strong'),family='gaussian')
inpartynewpid.dem.raw.cohort <- APCI::cohortdeviation(inpartynewpid.dem.raw$A,inpartynewpid.dem.raw$P,inpartynewpid.dem.raw$C,model=inpartynewpid.dem.raw$model,covariate=c('strong'))
## Model 3
inpartynewsort.dem.raw <- APCI::temp_model(data=anes.cumulative.dem,outcome='inparty.new',age='agegroup',period='yeargroup',cohort='cohort',covariate=c('ideosort','socialsort'),family='gaussian')
inpartynewsort.dem.raw.cohort <- APCI::cohortdeviation(inpartynewsort.dem.raw$A,inpartynewsort.dem.raw$P,inpartynewsort.dem.raw$C,model=inpartynewsort.dem.raw$model,covariate=c('strong'))
## Model 4
inpartynew.dem.controls <- APCI::temp_model(data=anes.cumulative.dem,outcome='inparty.new',age='agegroup',period='yeargroup',cohort='cohort',covariate=c('gender','black','education','ideo','attendance','south'),family='gaussian')
inpartynew.dem.controls.cohort <- APCI::cohortdeviation(inpartynew.dem.controls$A,inpartynew.dem.controls$P,inpartynew.dem.controls$C,model=inpartynew.dem.controls$model,covariate=c('gender','black','education','ideo','attendance','south'))
## Model 5
inpartynewpid.dem.controls <- APCI::temp_model(data=anes.cumulative.dem,outcome='inparty.new',age='agegroup',period='yeargroup',cohort='cohort',covariate=c('strong','gender','black','education','ideo','attendance','south'),family='gaussian')
inpartynewpid.dem.controls.cohort <- APCI::cohortdeviation(inpartynewpid.dem.controls$A,inpartynewpid.dem.controls$P,inpartynewpid.dem.controls$C,model=inpartynewpid.dem.controls$model,covariate=c('strong','gender','black','education','ideo','attendance','south'))
## Model 6
inpartynewsort.dem.controls <- APCI::temp_model(data=anes.cumulative.dem,outcome='inparty.new',age='agegroup',period='yeargroup',cohort='cohort',covariate=c('ideosort','socialsort','gender','black','education','ideo','attendance','south'),family='gaussian')
inpartynewsort.dem.controls.cohort <- APCI::cohortdeviation(inpartynewsort.dem.controls$A,inpartynewsort.dem.controls$P,inpartynewsort.dem.controls$C,model=inpartynewsort.dem.controls$model,covariate=c('ideosort','socialsort','gender','black','education','ideo','attendance','south'))

# Table A19
## Model 1
inpartynew.rep.raw <- APCI::temp_model(data=repnot90,outcome='inparty.new',age='agegroup',period='yeargroup',cohort='cohort',covariate=c(),family='gaussian')
inpartynew.rep.raw.cohort <- APCI::cohortdeviation(inpartynew.rep.raw$A,inpartynew.rep.raw$P,inpartynew.rep.raw$C,model=inpartynew.rep.raw$model,covariate=c())
## Model 2
inpartynewpid.rep.raw <- APCI::temp_model(data=repnot90,outcome='inparty.new',age='agegroup',period='yeargroup',cohort='cohort',covariate=c('strong'),family='gaussian')
inpartynewpid.rep.raw.cohort <- APCI::cohortdeviation(inpartynewpid.rep.raw$A,inpartynewpid.rep.raw$P,inpartynewpid.rep.raw$C,model=inpartynewpid.rep.raw$model,covariate=c('strong'))
## Model 3
inpartynewsort.rep.raw <- APCI::temp_model(data=repnot90,outcome='inparty.new',age='agegroup',period='yeargroup',cohort='cohort',covariate=c('ideosort','socialsort'),family='gaussian')
inpartynewsort.rep.raw.cohort <- APCI::cohortdeviation(inpartynewsort.rep.raw$A,inpartynewsort.rep.raw$P,inpartynewsort.rep.raw$C,model=inpartynewsort.rep.raw$model,covariate=c('strong'))
## Model 4
inpartynew.rep.controls <- APCI::temp_model(data=repnot90,outcome='inparty.new',age='agegroup',period='yeargroup',cohort='cohort',covariate=c('gender','black','education','ideo','attendance','south'),family='gaussian')
inpartynew.rep.controls.cohort <- APCI::cohortdeviation(inpartynew.rep.controls$A,inpartynew.rep.controls$P,inpartynew.rep.controls$C,model=inpartynew.rep.controls$model,covariate=c('gender','black','education','ideo','attendance','south'))
## Model 5
inpartynewpid.rep.controls <- APCI::temp_model(data=repnot90,outcome='inparty.new',age='agegroup',period='yeargroup',cohort='cohort',covariate=c('strong','gender','black','education','ideo','attendance','south'),family='gaussian')
inpartynewpid.rep.controls.cohort <- APCI::cohortdeviation(inpartynewpid.rep.controls$A,inpartynewpid.rep.controls$P,inpartynewpid.rep.controls$C,model=inpartynewpid.rep.controls$model,covariate=c('strong','gender','black','education','ideo','attendance','south'))
## Model 6
inpartynewsort.rep.controls <- APCI::temp_model(data=repnot90,outcome='inparty.new',age='agegroup',period='yeargroup',cohort='cohort',covariate=c('ideosort','socialsort','gender','black','education','ideo','attendance','south'),family='gaussian')
inpartynewsort.rep.controls.cohort <- APCI::cohortdeviation(inpartynewsort.rep.controls$A,inpartynewsort.rep.controls$P,inpartynewsort.rep.controls$C,model=inpartynewsort.rep.controls$model,covariate=c('ideosort','socialsort','gender','black','education','ideo','attendance','south'))

# Figure A21
appendixageinpartynewie.all <- data.frame(age=c(20,25,30,35,40,45,50,55,60,65,70,75,80,85,90),
                                          effect=c(-3.518,-3.755,-2.868,-2.316,-1.164,-0.791,-0.708,0.696,1.325,1.705,2.571,1.752,2.041,0.941,4.087),
                                          se=c(0.58,0.514,0.457,0.416,0.405,0.397,0.391,0.402,0.436,0.475,0.551,0.664,0.767,1.048,1.809))
appendixageinpartynewie.all$lci <- appendixageinpartynewie.all$effect - (1.96*appendixageinpartynewie.all$se)
appendixageinpartynewie.all$uci <- appendixageinpartynewie.all$effect + (1.96*appendixageinpartynewie.all$se)
appendixageinpartynewie.all$color <- ifelse(appendixageinpartynewie.all$lci<0 & appendixageinpartynewie.all$uci>0,"Gray50","Black")
appendixageinpartynewie.all$group <- 1

appendixageinpartynewapci.all <- data.frame(age=c(25,30,35,40,45,50,55,60,65,70,75,80,85,90),
                                            effect=c(summary(inpartynew.all.raw$model)$coefficients[2:15,1]),
                                            se=c(summary(inpartynew.all.raw$model)$coefficients[2:15,2]))
appendixageinpartynewapci.all$lci <- appendixageinpartynewapci.all$effect - (1.96*appendixageinpartynewapci.all$se)
appendixageinpartynewapci.all$uci <- appendixageinpartynewapci.all$effect + (1.96*appendixageinpartynewapci.all$se)
appendixageinpartynewapci.all$color <- ifelse(appendixageinpartynewapci.all$lci<0 & appendixageinpartynewapci.all$uci>0,"Gray50","Black")
appendixageinpartynewapci.all$group <- 1

figurea21.a <- ggplot(appendixageinpartynewapci.all,aes(x=age,y=effect,color=color,group=group)) + geom_point() + geom_errorbar(aes(ymin=lci,ymax=uci),size=0.5,width=0) + theme_classic() + theme(plot.title = element_text(hjust = 0.5)) + ggtitle("APCI") + geom_hline(yintercept=0,linetype="dashed",color="red") + scale_color_identity() + xlab("Age") + ylab("Deviation from Grand Mean") + ylim(-19,37)
figurea21.b <- ggplot(appendixageinpartynewie.all,aes(x=age,y=effect,color=color,group=group)) + geom_point() + geom_errorbar(aes(ymin=lci,ymax=uci),size=0.5,width=0) + theme_classic() + theme(plot.title = element_text(hjust = 0.5)) + ggtitle("IE") + geom_hline(yintercept=0,linetype="dashed",color="red") + scale_color_identity() + xlab("Age") + ylab("Deviation from Grand Mean") + ylim(-19,37)
ggarrange(figurea21.a,figurea21.b,nrow=1)

# Figure A22
appendixperiodinpartynewie.all <- data.frame(period=c(1975,1980,1985,1990,1995,2000,2005,2010,2015),
                                             effect=c(0.156,3.003,2.499,-1.331,0.073,1.582,1.071,-1.349,-5.702),
                                             se=c(0.404,0.377,0.294,0.301,0.275,0.507,0.403,0.337,0.421))
appendixperiodinpartynewie.all$lci <- appendixperiodinpartynewie.all$effect - (1.96*appendixperiodinpartynewie.all$se)
appendixperiodinpartynewie.all$uci <- appendixperiodinpartynewie.all$effect + (1.96*appendixperiodinpartynewie.all$se)
appendixperiodinpartynewie.all$color <- ifelse(appendixperiodinpartynewie.all$lci<0 & appendixperiodinpartynewie.all$uci>0,"Gray50","Black")
appendixperiodinpartynewie.all$group <- 1

appendixperiodinpartynewapci.all <- data.frame(period=c(1980,1985,1990,1995,2000,2005,2010,2015),
                                               effect=c(summary(inpartynew.all.raw$model)$coefficients[16:23,1]),
                                               se=c(summary(inpartynew.all.raw$model)$coefficients[16:23,2]))
appendixperiodinpartynewapci.all$lci <- appendixperiodinpartynewapci.all$effect - (1.96*appendixperiodinpartynewapci.all$se)
appendixperiodinpartynewapci.all$uci <- appendixperiodinpartynewapci.all$effect + (1.96*appendixperiodinpartynewapci.all$se)
appendixperiodinpartynewapci.all$color <- ifelse(appendixperiodinpartynewapci.all$lci<0 & appendixperiodinpartynewapci.all$uci>0,"Gray50","Black")
appendixperiodinpartynewapci.all$group <- 1

figurea22.a <- ggplot(appendixperiodinpartynewapci.all,aes(x=period,y=effect,color=color,group=group)) + geom_point() + geom_errorbar(aes(ymin=lci,ymax=uci),size=0.5,width=0) + theme_classic() + theme(plot.title = element_text(hjust = 0.5)) + ggtitle("APCI") + geom_hline(yintercept=0,linetype="dashed",color="red") + scale_color_identity() + xlab("Year") + ylab("Deviation from Grand Mean") + ylim(-19,37)
figurea22.b <- ggplot(appendixperiodinpartynewie.all,aes(x=period,y=effect,color=color,group=group)) + geom_point() + geom_errorbar(aes(ymin=lci,ymax=uci),size=0.5,width=0) + theme_classic() + theme(plot.title = element_text(hjust = 0.5)) + ggtitle("IE") + geom_hline(yintercept=0,linetype="dashed",color="red") + scale_color_identity() + xlab("Year") + ylab("Deviation from Grand Mean") + ylim(-19,37)
ggarrange(figurea22.a,figurea22.b,nrow=1)

# Figure A23
appendixcohortinpartynewie.all <- data.frame(cohort=c(1885,1890,1895,1900,1905,1910,1915,1920,1925,1930,1935,1940,1945,1950,1955,1960,1965,1970,1975,1980,1985,1990,1995),
                                             effect=c(0.642,10.097,0.99,0.908,0.918,-1.038,0.329,-0.840,-1.065,-1.112,-1.613,-2.162,-2.637,-2.138,-1.745,-1.230,-1.109,-0.697,-0.410,-0.215,0.41,1.790,1.928),
                                             se=c(5.717,3.27,1.921,1.519,1.179,1.026,0.928,0.841,0.753,0.707,0.634,0.551,0.471,0.406,0.354,0.345,0.384,0.459,0.528,0.608,0.660,0.759,1.228))
appendixcohortinpartynewie.all$lci <- appendixcohortinpartynewie.all$effect - (1.96*appendixcohortinpartynewie.all$se)
appendixcohortinpartynewie.all$uci <- appendixcohortinpartynewie.all$effect + (1.96*appendixcohortinpartynewie.all$se)
appendixcohortinpartynewie.all$color <- ifelse(appendixcohortinpartynewie.all$lci<0 & appendixcohortinpartynewie.all$uci>0,"Gray50","Black")
appendixcohortinpartynewie.all$group <- 1

appendixcohortinpartynewapci.all <- data.frame(cohort=c(1885,1890,1895,1900,1905,1910,1915,1920,1925,1930,1935,1940,1945,1950,1955,1960,1965,1970,1975,1980,1985,1990,1995),
                                               effect=c(as.numeric(inpartynew.all.raw.cohort$cohort_average$cohort_average)),
                                               se=c(as.numeric(inpartynew.all.raw.cohort$cohort_average$cohort_average_se)))
appendixcohortinpartynewapci.all$lci <- appendixcohortinpartynewapci.all$effect - (1.96*appendixcohortinpartynewapci.all$se)
appendixcohortinpartynewapci.all$uci <- appendixcohortinpartynewapci.all$effect + (1.96*appendixcohortinpartynewapci.all$se)
appendixcohortinpartynewapci.all$color <- ifelse(appendixcohortinpartynewapci.all$lci<0 & appendixcohortinpartynewapci.all$uci>0,"Gray50","Black")
appendixcohortinpartynewapci.all$group <- 1

figurea23.a <- ggplot(appendixcohortinpartynewapci.all,aes(x=cohort,y=effect,color=color,group=group)) + geom_point() + geom_errorbar(aes(ymin=lci,ymax=uci),size=0.5,width=0) + theme_classic() + theme(plot.title = element_text(hjust = 0.5)) + ggtitle("APCI") + geom_hline(yintercept=0,linetype="dashed",color="red") + scale_color_identity() + xlab("Birth Cohort") + ylab("Deviation from Grand Mean") + ylim(-19,37)
figurea23.b <- ggplot(appendixcohortinpartynewie.all,aes(x=cohort,y=effect,color=color,group=group)) + geom_point() + geom_errorbar(aes(ymin=lci,ymax=uci),size=0.5,width=0) + theme_classic() + theme(plot.title = element_text(hjust = 0.5)) + ggtitle("IE") + geom_hline(yintercept=0,linetype="dashed",color="red") + scale_color_identity() + xlab("Birth Cohort") + ylab("Deviation from Grand Mean") + ylim(-19,37)
ggarrange(figurea23.a,figurea23.b,nrow=1)

# Figure A24
appendixcohortinpartynewapcislopes.all <- data.frame(cohort=c(1885,1890,1895,1900,1905,1910,1915,1920,1925,1930,1935,1940,1945,1950,1955,1960,1965,1970,1975,1980,1985,1990,1995),
                                                     effect=c(as.numeric(inpartynew.all.raw.cohort$cohort_slope$cohort_slope)),
                                                     se=c(as.numeric(inpartynew.all.raw.cohort$cohort_slope$cohort_slope_se)))
appendixcohortinpartynewapcislopes.all$lci <- appendixcohortinpartynewapcislopes.all$effect - (1.96*appendixcohortinpartynewapcislopes.all$se)
appendixcohortinpartynewapcislopes.all$uci <- appendixcohortinpartynewapcislopes.all$effect + (1.96*appendixcohortinpartynewapcislopes.all$se)
appendixcohortinpartynewapcislopes.all$color <- ifelse(appendixcohortinpartynewapcislopes.all$lci<0 & appendixcohortinpartynewapcislopes.all$uci>0,"Gray50","Black")
appendixcohortinpartynewapcislopes.all$group <- 1

ggplot(appendixcohortinpartynewapcislopes.all,aes(x=cohort,y=effect,color=color,group=group)) + geom_point() + geom_errorbar(aes(ymin=lci,ymax=uci),size=0.5,width=0) + theme_classic() + theme(plot.title = element_text(hjust = 0.5)) + geom_hline(yintercept=0,linetype="dashed",color="red") + scale_color_identity() + xlab("Birth Cohort") + ylab("Deviation from Grand Mean") + ylim(-19,37)

# Figure A25
appendixageinpartynewie.dem <- data.frame(age=c(20,25,30,35,40,45,50,55,60,65,70,75,80,85,90),
                                          effect=c(-4.918,-5.069,-3.149,-2.684,-2.073,-0.167,-0.087,1.402,1.459,2.036,3.537,2.158,2.003,1.226,4.326),
                                          se=c(0.755,0.667,0.597,0.542,0.528,0.521,0.514,0.53,0.571,0.62,0.725,0.846,1.008,1.395,2.554))
appendixageinpartynewie.dem$lci <- appendixageinpartynewie.dem$effect - (1.96*appendixageinpartynewie.dem$se)
appendixageinpartynewie.dem$uci <- appendixageinpartynewie.dem$effect + (1.96*appendixageinpartynewie.dem$se)
appendixageinpartynewie.dem$color <- ifelse(appendixageinpartynewie.dem$lci<0 & appendixageinpartynewie.dem$uci>0,"Gray50","Black")
appendixageinpartynewie.dem$group <- 1

appendixageinpartynewapci.dem <- data.frame(age=c(25,30,35,40,45,50,55,60,65,70,75,80,85,90),
                                            effect=c(summary(inpartynew.dem.raw$model)$coefficients[2:15,1]),
                                            se=c(summary(inpartynew.dem.raw$model)$coefficients[2:15,2]))
appendixageinpartynewapci.dem$lci <- appendixageinpartynewapci.dem$effect - (1.96*appendixageinpartynewapci.dem$se)
appendixageinpartynewapci.dem$uci <- appendixageinpartynewapci.dem$effect + (1.96*appendixageinpartynewapci.dem$se)
appendixageinpartynewapci.dem$color <- ifelse(appendixageinpartynewapci.dem$lci<0 & appendixageinpartynewapci.dem$uci>0,"Gray50","Black")
appendixageinpartynewapci.dem$group <- 1

figurea25.a <- ggplot(appendixageinpartynewapci.dem,aes(x=age,y=effect,color=color,group=group)) + geom_point() + geom_errorbar(aes(ymin=lci,ymax=uci),size=0.5,width=0) + theme_classic() + theme(plot.title = element_text(hjust = 0.5)) + ggtitle("APCI") + geom_hline(yintercept=0,linetype="dashed",color="red") + scale_color_identity() + xlab("Age") + ylab("Deviation from Grand Mean") + ylim(-19,37)
figurea25.b <- ggplot(appendixageinpartynewie.dem,aes(x=age,y=effect,color=color,group=group)) + geom_point() + geom_errorbar(aes(ymin=lci,ymax=uci),size=0.5,width=0) + theme_classic() + theme(plot.title = element_text(hjust = 0.5)) + ggtitle("IE") + geom_hline(yintercept=0,linetype="dashed",color="red") + scale_color_identity() + xlab("Age") + ylab("Deviation from Grand Mean") + ylim(-19,37)
ggarrange(figurea25.a,figurea25.b,nrow=1)

# Figure A26
appendixperiodinpartynewie.dem <- data.frame(period=c(1975,1980,1985,1990,1995,2000,2005,2010,2015),
                                             effect=c(0.38,2.592,1.853,-1.816,0.003,-0.332,1.909,-0.628,-3.961),
                                             se=c(0.519,0.492,0.384,0.398,0.358,0.673,0.494,0.437,0.565))
appendixperiodinpartynewie.dem$lci <- appendixperiodinpartynewie.dem$effect - (1.96*appendixperiodinpartynewie.dem$se)
appendixperiodinpartynewie.dem$uci <- appendixperiodinpartynewie.dem$effect + (1.96*appendixperiodinpartynewie.dem$se)
appendixperiodinpartynewie.dem$color <- ifelse(appendixperiodinpartynewie.dem$lci<0 & appendixperiodinpartynewie.dem$uci>0,"Gray50","Black")
appendixperiodinpartynewie.dem$group <- 1

appendixperiodinpartynewapci.dem <- data.frame(period=c(1980,1985,1990,1995,2000,2005,2010,2015),
                                               effect=c(summary(inpartynew.dem.raw$model)$coefficients[16:23,1]),
                                               se=c(summary(inpartynew.dem.raw$model)$coefficients[16:23,2]))
appendixperiodinpartynewapci.dem$lci <- appendixperiodinpartynewapci.dem$effect - (1.96*appendixperiodinpartynewapci.dem$se)
appendixperiodinpartynewapci.dem$uci <- appendixperiodinpartynewapci.dem$effect + (1.96*appendixperiodinpartynewapci.dem$se)
appendixperiodinpartynewapci.dem$color <- ifelse(appendixperiodinpartynewapci.dem$lci<0 & appendixperiodinpartynewapci.dem$uci>0,"Gray50","Black")
appendixperiodinpartynewapci.dem$group <- 1

figurea26.a <- ggplot(appendixperiodinpartynewapci.dem,aes(x=period,y=effect,color=color,group=group)) + geom_point() + geom_errorbar(aes(ymin=lci,ymax=uci),size=0.5,width=0) + theme_classic() + theme(plot.title = element_text(hjust = 0.5)) + ggtitle("APCI") + geom_hline(yintercept=0,linetype="dashed",color="red") + scale_color_identity() + xlab("Year") + ylab("Deviation from Grand Mean") + ylim(-19,37)
figurea26.b <- ggplot(appendixperiodinpartynewie.dem,aes(x=period,y=effect,color=color,group=group)) + geom_point() + geom_errorbar(aes(ymin=lci,ymax=uci),size=0.5,width=0) + theme_classic() + theme(plot.title = element_text(hjust = 0.5)) + ggtitle("IE") + geom_hline(yintercept=0,linetype="dashed",color="red") + scale_color_identity() + xlab("Year") + ylab("Deviation from Grand Mean") + ylim(-19,37)
ggarrange(figurea26.a,figurea26.b,nrow=1)

# Figure A27
appendixcohortinpartynewie.dem <- data.frame(cohort=c(1885,1890,1895,1900,1905,1910,1915,1920,1925,1930,1935,1940,1945,1950,1955,1960,1965,1970,1975,1980,1985,1990,1995),
                                             effect=c(11.680,7.825,-1.283,0.079,-0.465,-2.318,-0.335,-1.319,-1.206,-2.21,-2.878,-2.397,-3.574,-2.65,-2.06,-1.023,-0.404,-0.493,0.616,-0.525,1.931,1.773,1.237),
                                             se=c(7.235,4.915,2.713,2.012,1.525,1.335,1.202,1.085,0.986,0.931,0.838,0.727,0.616,0.525,0.461,0.455,0.505,0.597,0.665,0.749,0.838,0.960,1.587))
appendixcohortinpartynewie.dem$lci <- appendixcohortinpartynewie.dem$effect - (1.96*appendixcohortinpartynewie.dem$se)
appendixcohortinpartynewie.dem$uci <- appendixcohortinpartynewie.dem$effect + (1.96*appendixcohortinpartynewie.dem$se)
appendixcohortinpartynewie.dem$color <- ifelse(appendixcohortinpartynewie.dem$lci<0 & appendixcohortinpartynewie.dem$uci>0,"Gray50","Black")
appendixcohortinpartynewie.dem$group <- 1

appendixcohortinpartynewapci.dem <- data.frame(cohort=c(1885,1890,1895,1900,1905,1910,1915,1920,1925,1930,1935,1940,1945,1950,1955,1960,1965,1970,1975,1980,1985,1990,1995),
                                               effect=c(as.numeric(inpartynew.dem.raw.cohort$cohort_average$cohort_average)),
                                               se=c(as.numeric(inpartynew.dem.raw.cohort$cohort_average$cohort_average_se)))
appendixcohortinpartynewapci.dem$lci <- appendixcohortinpartynewapci.dem$effect - (1.96*appendixcohortinpartynewapci.dem$se)
appendixcohortinpartynewapci.dem$uci <- appendixcohortinpartynewapci.dem$effect + (1.96*appendixcohortinpartynewapci.dem$se)
appendixcohortinpartynewapci.dem$color <- ifelse(appendixcohortinpartynewapci.dem$lci<0 & appendixcohortinpartynewapci.dem$uci>0,"Gray50","Black")
appendixcohortinpartynewapci.dem$group <- 1

figurea27.a <- ggplot(appendixcohortinpartynewapci.dem,aes(x=cohort,y=effect,color=color,group=group)) + geom_point() + geom_errorbar(aes(ymin=lci,ymax=uci),size=0.5,width=0) + theme_classic() + theme(plot.title = element_text(hjust = 0.5)) + ggtitle("APCI") + geom_hline(yintercept=0,linetype="dashed",color="red") + scale_color_identity() + xlab("Birth Cohort") + ylab("Deviation from Grand Mean") + ylim(-19,37)
figurea27.b <- ggplot(appendixcohortinpartynewie.dem,aes(x=cohort,y=effect,color=color,group=group)) + geom_point() + geom_errorbar(aes(ymin=lci,ymax=uci),size=0.5,width=0) + theme_classic() + theme(plot.title = element_text(hjust = 0.5)) + ggtitle("IE") + geom_hline(yintercept=0,linetype="dashed",color="red") + scale_color_identity() + xlab("Birth Cohort") + ylab("Deviation from Grand Mean") + ylim(-19,37)
ggarrange(figurea27.a,figurea27.b,nrow=1)

# Figure A28
appendixcohortinpartynewapcislopes.dem <- data.frame(cohort=c(1885,1890,1895,1900,1905,1910,1915,1920,1925,1930,1935,1940,1945,1950,1955,1960,1965,1970,1975,1980,1985,1990,1995),
                                                     effect=c(as.numeric(inpartynew.dem.raw.cohort$cohort_slope$cohort_slope)),
                                                     se=c(as.numeric(inpartynew.dem.raw.cohort$cohort_slope$cohort_slope_se)))
appendixcohortinpartynewapcislopes.dem$lci <- appendixcohortinpartynewapcislopes.dem$effect - (1.96*appendixcohortinpartynewapcislopes.dem$se)
appendixcohortinpartynewapcislopes.dem$uci <- appendixcohortinpartynewapcislopes.dem$effect + (1.96*appendixcohortinpartynewapcislopes.dem$se)
appendixcohortinpartynewapcislopes.dem$color <- ifelse(appendixcohortinpartynewapcislopes.dem$lci<0 & appendixcohortinpartynewapcislopes.dem$uci>0,"Gray50","Black")
appendixcohortinpartynewapcislopes.dem$group <- 1

ggplot(appendixcohortinpartynewapcislopes.dem,aes(x=cohort,y=effect,color=color,group=group)) + geom_point() + geom_errorbar(aes(ymin=lci,ymax=uci),size=0.5,width=0) + theme_classic() + theme(plot.title = element_text(hjust = 0.5)) + geom_hline(yintercept=0,linetype="dashed",color="red") + scale_color_identity() + xlab("Birth Cohort") + ylab("Deviation from Grand Mean") + ylim(-19,37)

# Figure A29
appendixageinpartynewie.rep <- data.frame(age=c(20,25,30,35,40,45,50,55,60,65,70,75,80,85,90),
                                          effect=c(-1.190,-1.890,-2.379,-1.9,-0.126,-1.953,-1.842,-0.512,0.997,0.923,1.350,1.319,2.073,0.494,4.638),
                                          se=c(0.902,0.801,0.708,0.642,0.623,0.612,0.593,0.604,0.663,0.726,0.834,0.969,1.16,1.56,2.554))
appendixageinpartynewie.rep$lci <- appendixageinpartynewie.rep$effect - (1.96*appendixageinpartynewie.rep$se)
appendixageinpartynewie.rep$uci <- appendixageinpartynewie.rep$effect + (1.96*appendixageinpartynewie.rep$se)
appendixageinpartynewie.rep$color <- ifelse(appendixageinpartynewie.rep$lci<0 & appendixageinpartynewie.rep$uci>0,"Gray50","Black")
appendixageinpartynewie.rep$group <- 1

appendixageinpartynewapci.rep <- data.frame(age=c(25,30,35,40,45,50,55,60,65,70,75,80,85),
                                            effect=c(summary(inpartynew.rep.raw$model)$coefficients[2:14,1]),
                                            se=c(summary(inpartynew.rep.raw$model)$coefficients[2:14,2]))
appendixageinpartynewapci.rep$lci <- appendixageinpartynewapci.rep$effect - (1.96*appendixageinpartynewapci.rep$se)
appendixageinpartynewapci.rep$uci <- appendixageinpartynewapci.rep$effect + (1.96*appendixageinpartynewapci.rep$se)
appendixageinpartynewapci.rep$color <- ifelse(appendixageinpartynewapci.rep$lci<0 & appendixageinpartynewapci.rep$uci>0,"Gray50","Black")
appendixageinpartynewapci.rep$group <- 1

figurea29.a <- ggplot(appendixageinpartynewapci.rep,aes(x=age,y=effect,color=color,group=group)) + geom_point() + geom_errorbar(aes(ymin=lci,ymax=uci),size=0.5,width=0) + theme_classic() + theme(plot.title = element_text(hjust = 0.5)) + ggtitle("APCI") + geom_hline(yintercept=0,linetype="dashed",color="red") + scale_color_identity() + xlab("Age") + ylab("Deviation from Grand Mean") + ylim(-19,37)
figurea29.b <- ggplot(appendixageinpartynewie.rep,aes(x=age,y=effect,color=color,group=group)) + geom_point() + geom_errorbar(aes(ymin=lci,ymax=uci),size=0.5,width=0) + theme_classic() + theme(plot.title = element_text(hjust = 0.5)) + ggtitle("IE") + geom_hline(yintercept=0,linetype="dashed",color="red") + scale_color_identity() + xlab("Age") + ylab("Deviation from Grand Mean") + ylim(-19,37)
ggarrange(figurea29.a,figurea29.b,nrow=1)

# Figure A30
appendixperiodinpartynewie.rep <- data.frame(period=c(1975,1980,1985,1990,1995,2000,2005,2010,2015),
                                             effect=c(-0.577,3.525,3.4,-0.44,0.357,4.538,-1.87,-2.461,-6.471),
                                             se=c(0.638,0.58,0.45,0.451,0.421,0.753,0.687,0.523,0.629))
appendixperiodinpartynewie.rep$lci <- appendixperiodinpartynewie.rep$effect - (1.96*appendixperiodinpartynewie.rep$se)
appendixperiodinpartynewie.rep$uci <- appendixperiodinpartynewie.rep$effect + (1.96*appendixperiodinpartynewie.rep$se)
appendixperiodinpartynewie.rep$color <- ifelse(appendixperiodinpartynewie.rep$lci<0 & appendixperiodinpartynewie.rep$uci>0,"Gray50","Black")
appendixperiodinpartynewie.rep$group <- 1

appendixperiodinpartynewapci.rep <- data.frame(period=c(1980,1985,1990,1995,2000,2005,2010,2015),
                                               effect=c(summary(inpartynew.rep.raw$model)$coefficients[15:22,1]),
                                               se=c(summary(inpartynew.rep.raw$model)$coefficients[15:22,2]))
appendixperiodinpartynewapci.rep$lci <- appendixperiodinpartynewapci.rep$effect - (1.96*appendixperiodinpartynewapci.rep$se)
appendixperiodinpartynewapci.rep$uci <- appendixperiodinpartynewapci.rep$effect + (1.96*appendixperiodinpartynewapci.rep$se)
appendixperiodinpartynewapci.rep$color <- ifelse(appendixperiodinpartynewapci.rep$lci<0 & appendixperiodinpartynewapci.rep$uci>0,"Gray50","Black")
appendixperiodinpartynewapci.rep$group <- 1

figurea30.a <- ggplot(appendixperiodinpartynewapci.rep,aes(x=period,y=effect,color=color,group=group)) + geom_point() + geom_errorbar(aes(ymin=lci,ymax=uci),size=0.5,width=0) + theme_classic() + theme(plot.title = element_text(hjust = 0.5)) + ggtitle("APCI") + geom_hline(yintercept=0,linetype="dashed",color="red") + scale_color_identity() + xlab("Year") + ylab("Deviation from Grand Mean") + ylim(-19,37)
figurea30.b <- ggplot(appendixperiodinpartynewie.rep,aes(x=period,y=effect,color=color,group=group)) + geom_point() + geom_errorbar(aes(ymin=lci,ymax=uci),size=0.5,width=0) + theme_classic() + theme(plot.title = element_text(hjust = 0.5)) + ggtitle("IE") + geom_hline(yintercept=0,linetype="dashed",color="red") + scale_color_identity() + xlab("Year") + ylab("Deviation from Grand Mean") + ylim(-19,37)
ggarrange(figurea30.a,figurea30.b,nrow=1)

# Figure A31
appendixcohortinpartynewie.rep <- data.frame(cohort=c(1885,1890,1895,1900,1905,1910,1915,1920,1925,1930,1935,1940,1945,1950,1955,1960,1965,1970,1975,1980,1985,1990,1995),
                                             effect=c(18.499,13.613,4.787,3.095,2.972,1.193,1.328,-0.211,-0.347,1.223,0.884,-0.961,-0.687,-1.087,-0.730,-0.982,-1.760,-0.405,-2.456,-0.795,-2.223,0.830,1.217),
                                             se=c(9.164,4.382,2.719,2.284,1.840,1.587,1.449,1.323,1.157,1.076,0.960,0.835,0.722,0.629,0.541,0.518,0.580,0.7,0.849,1.022,1.05,1.217,1.897))
appendixcohortinpartynewie.rep$lci <- appendixcohortinpartynewie.rep$effect - (1.96*appendixcohortinpartynewie.rep$se)
appendixcohortinpartynewie.rep$uci <- appendixcohortinpartynewie.rep$effect + (1.96*appendixcohortinpartynewie.rep$se)
appendixcohortinpartynewie.rep$color <- ifelse(appendixcohortinpartynewie.rep$lci<0 & appendixcohortinpartynewie.rep$uci>0,"Gray50","Black")
appendixcohortinpartynewie.rep$group <- 1

appendixcohortinpartynewapci.rep <- data.frame(cohort=c(1890,1895,1900,1905,1910,1915,1920,1925,1930,1935,1940,1945,1950,1955,1960,1965,1970,1975,1980,1985,1990,1995),
                                               effect=c(as.numeric(inpartynew.rep.raw.cohort$cohort_average$cohort_average)),
                                               se=c(as.numeric(inpartynew.rep.raw.cohort$cohort_average$cohort_average_se)))
appendixcohortinpartynewapci.rep$lci <- appendixcohortinpartynewapci.rep$effect - (1.96*appendixcohortinpartynewapci.rep$se)
appendixcohortinpartynewapci.rep$uci <- appendixcohortinpartynewapci.rep$effect + (1.96*appendixcohortinpartynewapci.rep$se)
appendixcohortinpartynewapci.rep$color <- ifelse(appendixcohortinpartynewapci.rep$lci<0 & appendixcohortinpartynewapci.rep$uci>0,"Gray50","Black")
appendixcohortinpartynewapci.rep$group <- 1

figurea31.a <- ggplot(appendixcohortinpartynewapci.rep,aes(x=cohort,y=effect,color=color,group=group)) + geom_point() + geom_errorbar(aes(ymin=lci,ymax=uci),size=0.5,width=0) + theme_classic() + theme(plot.title = element_text(hjust = 0.5)) + ggtitle("APCI") + geom_hline(yintercept=0,linetype="dashed",color="red") + scale_color_identity() + xlab("Birth Cohort") + ylab("Deviation from Grand Mean") + ylim(-19,37)
figurea31.b <- ggplot(appendixcohortinpartynewie.rep,aes(x=cohort,y=effect,color=color,group=group)) + geom_point() + geom_errorbar(aes(ymin=lci,ymax=uci),size=0.5,width=0) + theme_classic() + theme(plot.title = element_text(hjust = 0.5)) + ggtitle("IE") + geom_hline(yintercept=0,linetype="dashed",color="red") + scale_color_identity() + xlab("Birth Cohort") + ylab("Deviation from Grand Mean") + ylim(-19,37)
ggarrange(figurea31.a,figurea31.b,nrow=1)

# Figure A32
appendixcohortinpartynewapcislopes.rep <- data.frame(cohort=c(1890,1895,1900,1905,1910,1915,1920,1925,1930,1935,1940,1945,1950,1955,1960,1965,1970,1975,1980,1985,1990,1995),
                                                     effect=c(as.numeric(inpartynew.rep.raw.cohort$cohort_slope$cohort_slope)),
                                                     se=c(as.numeric(inpartynew.rep.raw.cohort$cohort_slope$cohort_slope_se)))
appendixcohortinpartynewapcislopes.rep$lci <- appendixcohortinpartynewapcislopes.rep$effect - (1.96*appendixcohortinpartynewapcislopes.rep$se)
appendixcohortinpartynewapcislopes.rep$uci <- appendixcohortinpartynewapcislopes.rep$effect + (1.96*appendixcohortinpartynewapcislopes.rep$se)
appendixcohortinpartynewapcislopes.rep$color <- ifelse(appendixcohortinpartynewapcislopes.rep$lci<0 & appendixcohortinpartynewapcislopes.rep$uci>0,"Gray50","Black")
appendixcohortinpartynewapcislopes.rep$group <- 1

ggplot(appendixcohortinpartynewapcislopes.rep,aes(x=cohort,y=effect,color=color,group=group)) + geom_point() + geom_errorbar(aes(ymin=lci,ymax=uci),size=0.5,width=0) + theme_classic() + theme(plot.title = element_text(hjust = 0.5)) + geom_hline(yintercept=0,linetype="dashed",color="red") + scale_color_identity() + xlab("Birth Cohort") + ylab("Deviation from Grand Mean") + ylim(-19,37)

# Table A21
allnot90 <- subset(anes.cumulative.full,agegroup!="90")
allnot90$agegroup <- droplevels(allnot90$agegroup)
## Model 1
inpartycombined.all.raw <- APCI::temp_model(data=anes.cumulative.full,outcome='inparty.combined',age='agegroup',period='yeargroup.old',cohort='cohort.combined',covariate=c(),family='gaussian')
inpartycombined.all.raw.cohort <- APCI::cohortdeviation(inpartycombined.all.raw$A,inpartycombined.all.raw$P,inpartycombined.all.raw$C,model=inpartycombined.all.raw$model,covariate=c())
## Model 2
inpartycombinedpid.all.raw <- APCI::temp_model(data=anes.cumulative.full,outcome='inparty.combined',age='agegroup',period='yeargroup.old',cohort='cohort.combined',covariate=c('strong'),family='gaussian')
inpartycombinedpid.all.raw.cohort <- APCI::cohortdeviation(inpartycombinedpid.all.raw$A,inpartycombinedpid.all.raw$P,inpartycombinedpid.all.raw$C,model=inpartycombinedpid.all.raw$model,covariate=c('strong'))
## Model 3
inpartycombinedsort.all.raw <- APCI::temp_model(data=anes.cumulative.full,outcome='inparty.combined',age='agegroup',period='yeargroup.sort',cohort='cohort.combined',covariate=c('ideosort','socialsort'),family='gaussian')
inpartycombinedsort.all.raw.cohort <- APCI::cohortdeviation(inpartycombinedsort.all.raw$A,inpartycombinedsort.all.raw$P,inpartycombinedsort.all.raw$C,model=inpartycombinedsort.all.raw$model,covariate=c('strong'))
## Model 4
inpartycombined.all.controls <- APCI::temp_model(data=anes.cumulative.full,outcome='inparty.combined',age='agegroup',period='yeargroup.sort',cohort='cohort.combined',covariate=c('gender','black','education','ideo','attendance','south'),family='gaussian')
inpartycombined.all.controls.cohort <- APCI::cohortdeviation(inpartycombined.all.controls$A,inpartycombined.all.controls$P,inpartycombined.all.controls$C,model=inpartycombined.all.controls$model,covariate=c('gender','black','education','ideo','attendance','south'))
## Model 5
inpartycombinedpid.all.controls <- APCI::temp_model(data=anes.cumulative.full,outcome='inparty.combined',age='agegroup',period='yeargroup.sort',cohort='cohort.combined',covariate=c('strong','gender','black','education','ideo','attendance','south'),family='gaussian')
inpartycombinedpid.all.controls.cohort <- APCI::cohortdeviation(inpartycombinedpid.all.controls$A,inpartycombinedpid.all.controls$P,inpartycombinedpid.all.controls$C,model=inpartycombinedpid.all.controls$model,covariate=c('strong','gender','black','education','ideo','attendance','south'))
## Model 6
inpartycombinedsort.all.controls <- APCI::temp_model(data=anes.cumulative.full,outcome='inparty.combined',age='agegroup',period='yeargroup.sort',cohort='cohort.combined',covariate=c('ideosort','socialsort','gender','black','education','ideo','attendance','south'),family='gaussian')
inpartycombinedsort.all.controls.cohort <- APCI::cohortdeviation(inpartycombinedsort.all.controls$A,inpartycombinedsort.all.controls$P,inpartycombinedsort.all.controls$C,model=inpartycombinedsort.all.controls$model,covariate=c('ideosort','socialsort','gender','black','education','ideo','attendance','south'))

# Table A22
demnot90 <- subset(anes.cumulative.dem,agegroup!="90")
demnot90$agegroup <- as.factor(ifelse(demnot90$agegroup=="20","20",ifelse(demnot90$agegroup=="25","25",ifelse(demnot90$agegroup=="30","30",ifelse(demnot90$agegroup=="35","35",ifelse(demnot90$agegroup=="40","40",ifelse(demnot90$agegroup=="45","45",ifelse(demnot90$agegroup=="50","50",ifelse(demnot90$agegroup=="55","55",ifelse(demnot90$agegroup=="60","60",ifelse(demnot90$agegroup=="65","65",ifelse(demnot90$agegroup=="70","70",ifelse(demnot90$agegroup=="75","75",ifelse(demnot90$agegroup=="80","80","85"))))))))))))))
demnot90$cohort.old <- as.factor(ifelse(demnot90$cohort.old=="1890","1890",ifelse(demnot90$cohort.old=="1895","1895",ifelse(demnot90$cohort.old=="1900","1900",ifelse(demnot90$cohort.old=="1905","1905",ifelse(demnot90$cohort.old=="1910","1910",ifelse(demnot90$cohort.old=="1915","1915",ifelse(demnot90$cohort.old=="1920","1920",ifelse(demnot90$cohort.old=="1925","1925",ifelse(demnot90$cohort.old=="1930","1930",ifelse(demnot90$cohort.old=="1935","1935",ifelse(demnot90$cohort.old=="1940","1940",ifelse(demnot90$cohort.old=="1945","1945",ifelse(demnot90$cohort.old=="1950","1950",ifelse(demnot90$cohort.old=="1955","1955",ifelse(demnot90$cohort.old=="1960","1960","1965"))))))))))))))))
demnot90$yeargroup.combined <- demnot90$yeargroup.old
demnot90$yeargroup.combined[demnot90$yeargroup.combined=="1965" | demnot90$yeargroup.combined=="1970"] <- NA
demnot90$yeargroup.combined <- droplevels(demnot90$yeargroup.combined)
demnot90$yeargroup.old <- as.factor(ifelse(demnot90$yeargroup.old=="1975","1975",ifelse(demnot90$yeargroup.old=="1980","1980",ifelse(demnot90$yeargroup.old=="1985","1985",NA))))
demnot90$yeargroup.sort <- as.factor(ifelse(demnot90$yeargroup.old=="1975","1975",ifelse(demnot90$yeargroup.old=="1980","1980",ifelse(demnot90$yeargroup.old=="1985","1985",NA))))
## Model 1
inpartycombined.dem.raw <- APCI::temp_model(data=anes.cumulative.dem,outcome='inparty.combined',age='agegroup',period='yeargroup.old',cohort='cohort.combined',covariate=c(),family='gaussian')
inpartycombined.dem.raw.cohort <- APCI::cohortdeviation(inpartycombined.dem.raw$A,inpartycombined.dem.raw$P,inpartycombined.dem.raw$C,model=inpartycombined.dem.raw$model,covariate=c())
## Model 2
inpartycombinedpid.dem.raw <- APCI::temp_model(data=anes.cumulative.dem,outcome='inparty.combined',age='agegroup',period='yeargroup.old',cohort='cohort.combined',covariate=c('strong'),family='gaussian')
inpartycombinedpid.dem.raw.cohort <- APCI::cohortdeviation(inpartycombinedpid.dem.raw$A,inpartycombinedpid.dem.raw$P,inpartycombinedpid.dem.raw$C,model=inpartycombinedpid.dem.raw$model,covariate=c('strong'))
## Model 3
inpartycombinedsort.dem.raw <- APCI::temp_model(data=demnot90,outcome='inparty.combined',age='agegroup',period='yeargroup.combined',cohort='cohort.combined',covariate=c('ideosort','socialsort'),family='gaussian')
inpartycombinedsort.dem.raw.cohort <- APCI::cohortdeviation(inpartycombinedsort.dem.raw$A,inpartycombinedsort.dem.raw$P,inpartycombinedsort.dem.raw$C,model=inpartycombinedsort.dem.raw$model,covariate=c('strong'))
## Model 4
inpartycombined.dem.controls <- APCI::temp_model(data=demnot90,outcome='inparty.combined',age='agegroup',period='yeargroup.combined',cohort='cohort.combined',covariate=c('gender','black','education','ideo','attendance','south'),family='gaussian')
inpartycombined.dem.controls.cohort <- APCI::cohortdeviation(inpartycombined.dem.controls$A,inpartycombined.dem.controls$P,inpartycombined.dem.controls$C,model=inpartycombined.dem.controls$model,covariate=c('gender','black','education','ideo','attendance','south'))
## Model 5
inpartycombinedpid.dem.controls <- APCI::temp_model(data=demnot90,outcome='inparty.combined',age='agegroup',period='yeargroup.combined',cohort='cohort.combined',covariate=c('strong','gender','black','education','ideo','attendance','south'),family='gaussian')
inpartycombinedpid.dem.controls.cohort <- APCI::cohortdeviation(inpartycombinedpid.dem.controls$A,inpartycombinedpid.dem.controls$P,inpartycombinedpid.dem.controls$C,model=inpartycombinedpid.dem.controls$model,covariate=c('strong','gender','black','education','ideo','attendance','south'))
## Model 6
inpartycombinedsort.dem.controls <- APCI::temp_model(data=demnot90,outcome='inparty.combined',age='agegroup',period='yeargroup.combined',cohort='cohort.combined',covariate=c('ideosort','socialsort','gender','black','education','ideo','attendance','south'),family='gaussian')
inpartycombinedsort.dem.controls.cohort <- APCI::cohortdeviation(inpartycombinedsort.dem.controls$A,inpartycombinedsort.dem.controls$P,inpartycombinedsort.dem.controls$C,model=inpartycombinedsort.dem.controls$model,covariate=c('ideosort','socialsort','gender','black','education','ideo','attendance','south'))

# Table A25
## Model 1
inpartycombined.rep.raw <- APCI::temp_model(data=repnot90,outcome='inparty.combined',age='agegroup',period='yeargroup.old',cohort='cohort.combined',covariate=c(),family='gaussian')
inpartycombined.rep.raw.cohort <- APCI::cohortdeviation(inpartycombined.rep.raw$A,inpartycombined.rep.raw$P,inpartycombined.rep.raw$C,model=inpartycombined.rep.raw$model,covariate=c())
## Model 2
inpartycombinedpid.rep.raw <- APCI::temp_model(data=repnot90,outcome='inparty.combined',age='agegroup',period='yeargroup.old',cohort='cohort.combined',covariate=c('strong'),family='gaussian')
inpartycombinedpid.rep.raw.cohort <- APCI::cohortdeviation(inpartycombinedpid.rep.raw$A,inpartycombinedpid.rep.raw$P,inpartycombinedpid.rep.raw$C,model=inpartycombinedpid.rep.raw$model,covariate=c('strong'))
## Model 3
inpartycombinedsort.rep.raw <- APCI::temp_model(data=repnot90,outcome='inparty.combined',age='agegroup',period='yeargroup.combined',cohort='cohort.combined',covariate=c('ideosort','socialsort'),family='gaussian')
inpartycombinedsort.rep.raw.cohort <- APCI::cohortdeviation(inpartycombinedsort.rep.raw$A,inpartycombinedsort.rep.raw$P,inpartycombinedsort.rep.raw$C,model=inpartycombinedsort.rep.raw$model,covariate=c('strong'))
## Model 4
inpartycombined.rep.controls <- APCI::temp_model(data=repnot90,outcome='inparty.combined',age='agegroup',period='yeargroup.combined',cohort='cohort.combined',covariate=c('gender','black','education','ideo','attendance','south'),family='gaussian')
inpartycombined.rep.controls.cohort <- APCI::cohortdeviation(inpartycombined.rep.controls$A,inpartycombined.rep.controls$P,inpartycombined.rep.controls$C,model=inpartycombined.rep.controls$model,covariate=c('gender','black','education','ideo','attendance','south'))
## Model 5
inpartycombinedpid.rep.controls <- APCI::temp_model(data=repnot90,outcome='inparty.combined',age='agegroup',period='yeargroup.combined',cohort='cohort.combined',covariate=c('strong','gender','black','education','ideo','attendance','south'),family='gaussian')
inpartycombinedpid.rep.controls.cohort <- APCI::cohortdeviation(inpartycombinedpid.rep.controls$A,inpartycombinedpid.rep.controls$P,inpartycombinedpid.rep.controls$C,model=inpartycombinedpid.rep.controls$model,covariate=c('strong','gender','black','education','ideo','attendance','south'))
## Model 6
inpartycombinedsort.rep.controls <- APCI::temp_model(data=repnot90,outcome='inparty.combined',age='agegroup',period='yeargroup.combined',cohort='cohort.combined',covariate=c('ideosort','socialsort','gender','black','education','ideo','attendance','south'),family='gaussian')
inpartycombinedsort.rep.controls.cohort <- APCI::cohortdeviation(inpartycombinedsort.rep.controls$A,inpartycombinedsort.rep.controls$P,inpartycombinedsort.rep.controls$C,model=inpartycombinedsort.rep.controls$model,covariate=c('ideosort','socialsort','gender','black','education','ideo','attendance','south'))

# Figure A33
appendixageinpartycombinedie.all <- data.frame(age=c(20,25,30,35,40,45,50,55,60,65,70,75,80,85,90),
                                               effect=c(-2.801,-2.825,-2.376,-1.790,-1.177,-0.811,-0.535,0.365,0.776,1.436,2.294,0.995,1.461,0.349,4.639),
                                               se=c(0.623,0.545,0.479,0.420,0.381,0.353,0.335,0.341,0.372,0.417,0.493,0.587,0.72,0.996,1.64))
appendixageinpartycombinedie.all$lci <- appendixageinpartycombinedie.all$effect - (1.96*appendixageinpartycombinedie.all$se)
appendixageinpartycombinedie.all$uci <- appendixageinpartycombinedie.all$effect + (1.96*appendixageinpartycombinedie.all$se)
appendixageinpartycombinedie.all$color <- ifelse(appendixageinpartycombinedie.all$lci<0 & appendixageinpartycombinedie.all$uci>0,"Gray50","Black")
appendixageinpartycombinedie.all$group <- 1

appendixageinpartycombinedapci.all <- data.frame(age=c(25,30,35,40,45,50,55,60,65,70,75,80,85,90),
                                                 effect=c(summary(inpartycombined.all.raw$model)$coefficients[2:15,1]),
                                                 se=c(summary(inpartycombined.all.raw$model)$coefficients[2:15,2]))
appendixageinpartycombinedapci.all$lci <- appendixageinpartycombinedapci.all$effect - (1.96*appendixageinpartycombinedapci.all$se)
appendixageinpartycombinedapci.all$uci <- appendixageinpartycombinedapci.all$effect + (1.96*appendixageinpartycombinedapci.all$se)
appendixageinpartycombinedapci.all$color <- ifelse(appendixageinpartycombinedapci.all$lci<0 & appendixageinpartycombinedapci.all$uci>0,"Gray50","Black")
appendixageinpartycombinedapci.all$group <- 1

figurea33.a <- ggplot(appendixageinpartycombinedapci.all,aes(x=age,y=effect,color=color,group=group)) + geom_point() + geom_errorbar(aes(ymin=lci,ymax=uci),size=0.5,width=0) + theme_classic() + theme(plot.title = element_text(hjust = 0.5)) + ggtitle("APCI") + geom_hline(yintercept=0,linetype="dashed",color="red") + scale_color_identity() + xlab("Age") + ylab("Deviation from Grand Mean") + ylim(-34,39)
figurea33.b <- ggplot(appendixageinpartycombinedie.all,aes(x=age,y=effect,color=color,group=group)) + geom_point() + geom_errorbar(aes(ymin=lci,ymax=uci),size=0.5,width=0) + theme_classic() + theme(plot.title = element_text(hjust = 0.5)) + ggtitle("IE") + geom_hline(yintercept=0,linetype="dashed",color="red") + scale_color_identity() + xlab("Age") + ylab("Deviation from Grand Mean") + ylim(-34,39)
ggarrange(figurea33.a,figurea33.b,nrow=1)

# Figure A34
appendixperiodinpartycombinedie.all <- data.frame(period=c(1960,1965,1970,1975,1980,1985,1990,1995,2000,2005,2010,2015),
                                                  effect=c(7.334,3.515,-0.438,-2.266,1.567,1.159,-2.555,-1.003,0.648,0.272,-1.997,-6.236),
                                                  se=c(0.601,0.464,0.426,0.338,0.338,0.256,0.292,0.298,0.545,0.473,0.446,0.543))
appendixperiodinpartycombinedie.all$lci <- appendixperiodinpartycombinedie.all$effect - (1.96*appendixperiodinpartycombinedie.all$se)
appendixperiodinpartycombinedie.all$uci <- appendixperiodinpartycombinedie.all$effect + (1.96*appendixperiodinpartycombinedie.all$se)
appendixperiodinpartycombinedie.all$color <- ifelse(appendixperiodinpartycombinedie.all$lci<0 & appendixperiodinpartycombinedie.all$uci>0,"Gray50","Black")
appendixperiodinpartycombinedie.all$group <- 1

appendixperiodinpartycombinedapci.all <- data.frame(period=c(1965,1970,1975,1980,1985,1990,1995,2000,2005,2010,2015),
                                                    effect=c(summary(inpartycombined.all.raw$model)$coefficients[16:26,1]),
                                                    se=c(summary(inpartycombined.all.raw$model)$coefficients[16:26,2]))
appendixperiodinpartycombinedapci.all$lci <- appendixperiodinpartycombinedapci.all$effect - (1.96*appendixperiodinpartycombinedapci.all$se)
appendixperiodinpartycombinedapci.all$uci <- appendixperiodinpartycombinedapci.all$effect + (1.96*appendixperiodinpartycombinedapci.all$se)
appendixperiodinpartycombinedapci.all$color <- ifelse(appendixperiodinpartycombinedapci.all$lci<0 & appendixperiodinpartycombinedapci.all$uci>0,"Gray50","Black")
appendixperiodinpartycombinedapci.all$group <- 1

figurea34.a <- ggplot(appendixperiodinpartycombinedapci.all,aes(x=period,y=effect,color=color,group=group)) + geom_point() + geom_errorbar(aes(ymin=lci,ymax=uci),size=0.5,width=0) + theme_classic() + theme(plot.title = element_text(hjust = 0.5)) + ggtitle("APCI") + geom_hline(yintercept=0,linetype="dashed",color="red") + scale_color_identity() + xlab("Year") + ylab("Deviation from Grand Mean") + ylim(-34,39)
figurea34.b <- ggplot(appendixperiodinpartycombinedie.all,aes(x=period,y=effect,color=color,group=group)) + geom_point() + geom_errorbar(aes(ymin=lci,ymax=uci),size=0.5,width=0) + theme_classic() + theme(plot.title = element_text(hjust = 0.5)) + ggtitle("IE") + geom_hline(yintercept=0,linetype="dashed",color="red") + scale_color_identity() + xlab("Year") + ylab("Deviation from Grand Mean") + ylim(-34,39)
ggarrange(figurea34.a,figurea34.b,nrow=1)

# Table A35
appendixcohortinpartycombinedie.all <- data.frame(cohort=c(1870,1875,1880,1885,1890,1895,1900,1905,1910,1915,1920,1925,1930,1935,1940,1945,1950,1955,1960,1965,1970,1975,1980,1985,1990,1995),
                                                  effect=c(5.177,-9.131,3.362,2.465,3.167,2.55,1.341,1.725,0.745,1.146,0.032,-0.098,-0.174,-0.910,-1.596,-2.144,-1.958,-1.800,-1.344,-1.301,-0.977,-0.855,-0.882,-0.408,0.859,1.01),
                                                  se=c(7.875,5.371,2.998,1.845,1.512,1.289,1.168,1.042,0.939,0.852,0.769,0.686,0.623,0.551,0.471,0.4,0.348,0.32,0.329,0.385,0.472,0.556,0.646,0.712,0.818,1.261))
appendixcohortinpartycombinedie.all$lci <- appendixcohortinpartycombinedie.all$effect - (1.96*appendixcohortinpartycombinedie.all$se)
appendixcohortinpartycombinedie.all$uci <- appendixcohortinpartycombinedie.all$effect + (1.96*appendixcohortinpartycombinedie.all$se)
appendixcohortinpartycombinedie.all$color <- ifelse(appendixcohortinpartycombinedie.all$lci<0 & appendixcohortinpartycombinedie.all$uci>0,"Gray50","Black")
appendixcohortinpartycombinedie.all$group <- 1

appendixcohortinpartycombinedapci.all <- data.frame(cohort=c(1870,1875,1880,1885,1890,1895,1900,1905,1910,1915,1920,1925,1930,1935,1940,1945,1950,1955,1960,1965,1970,1975,1980,1985,1990,1995),
                                                    effect=c(as.numeric(inpartycombined.all.raw.cohort$cohort_average$cohort_average)),
                                                    se=c(as.numeric(inpartycombined.all.raw.cohort$cohort_average$cohort_average_se)))
appendixcohortinpartycombinedapci.all$lci <- appendixcohortinpartycombinedapci.all$effect - (1.96*appendixcohortinpartycombinedapci.all$se)
appendixcohortinpartycombinedapci.all$uci <- appendixcohortinpartycombinedapci.all$effect + (1.96*appendixcohortinpartycombinedapci.all$se)
appendixcohortinpartycombinedapci.all$color <- ifelse(appendixcohortinpartycombinedapci.all$lci<0 & appendixcohortinpartycombinedapci.all$uci>0,"Gray50","Black")
appendixcohortinpartycombinedapci.all$group <- 1

figurea35.a <- ggplot(appendixcohortinpartycombinedapci.all,aes(x=cohort,y=effect,color=color,group=group)) + geom_point() + geom_errorbar(aes(ymin=lci,ymax=uci),size=0.5,width=0) + theme_classic() + theme(plot.title = element_text(hjust = 0.5)) + ggtitle("APCI") + geom_hline(yintercept=0,linetype="dashed",color="red") + scale_color_identity() + xlab("Birth Cohort") + ylab("Deviation from Grand Mean") + ylim(-34,39)
figurea35.b <- ggplot(appendixcohortinpartycombinedie.all,aes(x=cohort,y=effect,color=color,group=group)) + geom_point() + geom_errorbar(aes(ymin=lci,ymax=uci),size=0.5,width=0) + theme_classic() + theme(plot.title = element_text(hjust = 0.5)) + ggtitle("IE") + geom_hline(yintercept=0,linetype="dashed",color="red") + scale_color_identity() + xlab("Birth Cohort") + ylab("Deviation from Grand Mean") + ylim(-34,39)
ggarrange(figurea35.a,figurea35.b,nrow=1)

# Figure A36
appendixcohortinpartycombinedapcislopes.all <- data.frame(cohort=c(1870,1875,1880,1885,1890,1895,1900,1905,1910,1915,1920,1925,1930,1935,1940,1945,1950,1955,1960,1965,1970,1975,1980,1985,1990,1995),
                                                          effect=c(as.numeric(inpartycombined.all.raw.cohort$cohort_slope$cohort_slope)),
                                                          se=c(as.numeric(inpartycombined.all.raw.cohort$cohort_slope$cohort_slope_se)))
appendixcohortinpartycombinedapcislopes.all$lci <- appendixcohortinpartycombinedapcislopes.all$effect - (1.96*appendixcohortinpartycombinedapcislopes.all$se)
appendixcohortinpartycombinedapcislopes.all$uci <- appendixcohortinpartycombinedapcislopes.all$effect + (1.96*appendixcohortinpartycombinedapcislopes.all$se)
appendixcohortinpartycombinedapcislopes.all$color <- ifelse(appendixcohortinpartycombinedapcislopes.all$lci<0 & appendixcohortinpartycombinedapcislopes.all$uci>0,"Gray50","Black")
appendixcohortinpartycombinedapcislopes.all$group <- 1

ggplot(appendixcohortinpartycombinedapcislopes.all,aes(x=cohort,y=effect,color=color,group=group)) + geom_point() + geom_errorbar(aes(ymin=lci,ymax=uci),size=0.5,width=0) + theme_classic() + theme(plot.title = element_text(hjust = 0.5)) + geom_hline(yintercept=0,linetype="dashed",color="red") + scale_color_identity() + xlab("Birth Cohort") + ylab("Deviation from Grand Mean") + ylim(-34,39)

# Figure A37
appendixageinpartycombinedie.dem <- data.frame(age=c(20,25,30,35,40,45,50,55,60,65,70,75,80,85,90),
                                               effect=c(-3.429,-3.477,-2.405,-1.815,-1.747,-0.659,-0.267,0.594,0.742,1.38,3.087,1.001,1.202,0.365,5.428),
                                               se=c(1.069,0.929,0.8,0.679,0.581,0.502,0.449,0.446,0.498,0.580,0.715,0.876,1.078,1.455,2.428))
appendixageinpartycombinedie.dem$lci <- appendixageinpartycombinedie.dem$effect - (1.96*appendixageinpartycombinedie.dem$se)
appendixageinpartycombinedie.dem$uci <- appendixageinpartycombinedie.dem$effect + (1.96*appendixageinpartycombinedie.dem$se)
appendixageinpartycombinedie.dem$color <- ifelse(appendixageinpartycombinedie.dem$lci<0 & appendixageinpartycombinedie.dem$uci>0,"Gray50","Black")
appendixageinpartycombinedie.dem$group <- 1

appendixageinpartycombinedapci.dem <- data.frame(age=c(25,30,35,40,45,50,55,60,65,70,75,80,85,90),
                                                 effect=c(summary(inpartycombined.dem.raw$model)$coefficients[2:15,1]),
                                                 se=c(summary(inpartycombined.dem.raw$model)$coefficients[2:15,2]))
appendixageinpartycombinedapci.dem$lci <- appendixageinpartycombinedapci.dem$effect - (1.96*appendixageinpartycombinedapci.dem$se)
appendixageinpartycombinedapci.dem$uci <- appendixageinpartycombinedapci.dem$effect + (1.96*appendixageinpartycombinedapci.dem$se)
appendixageinpartycombinedapci.dem$color <- ifelse(appendixageinpartycombinedapci.dem$lci<0 & appendixageinpartycombinedapci.dem$uci>0,"Gray50","Black")
appendixageinpartycombinedapci.dem$group <- 1

figurea37.a <- ggplot(appendixageinpartycombinedapci.dem,aes(x=age,y=effect,color=color,group=group)) + geom_point() + geom_errorbar(aes(ymin=lci,ymax=uci),size=0.5,width=0) + theme_classic() + theme(plot.title = element_text(hjust = 0.5)) + ggtitle("APCI") + geom_hline(yintercept=0,linetype="dashed",color="red") + scale_color_identity() + xlab("Age") + ylab("Deviation from Grand Mean") + ylim(-34,39)
figurea37.b <- ggplot(appendixageinpartycombinedie.dem,aes(x=age,y=effect,color=color,group=group)) + geom_point() + geom_errorbar(aes(ymin=lci,ymax=uci),size=0.5,width=0) + theme_classic() + theme(plot.title = element_text(hjust = 0.5)) + ggtitle("IE") + geom_hline(yintercept=0,linetype="dashed",color="red") + scale_color_identity() + xlab("Age") + ylab("Deviation from Grand Mean") + ylim(-34,39)
ggarrange(figurea37.a,figurea37.b,nrow=1)

# Figure A38
appendixperiodinpartycombinedie.dem <- data.frame(period=c(1960,1965,1970,1975,1980,1985,1990,1995,2000,2005,2010,2015),
                                                  effect=c(7.458,2.787,-0.292,-2.218,0.892,0.324,-3.134,-1.064,-1.162,1.349,-0.913,-4.026),
                                                  se=c(0.939,0.748,0.653,0.5,0.464,0.338,0.391,0.416,0.768,0.687,0.735,0.907))
appendixperiodinpartycombinedie.dem$lci <- appendixperiodinpartycombinedie.dem$effect - (1.96*appendixperiodinpartycombinedie.dem$se)
appendixperiodinpartycombinedie.dem$uci <- appendixperiodinpartycombinedie.dem$effect + (1.96*appendixperiodinpartycombinedie.dem$se)
appendixperiodinpartycombinedie.dem$color <- ifelse(appendixperiodinpartycombinedie.dem$lci<0 & appendixperiodinpartycombinedie.dem$uci>0,"Gray50","Black")
appendixperiodinpartycombinedie.dem$group <- 1

appendixperiodinpartycombinedapci.dem <- data.frame(period=c(1965,1970,1975,1980,1985,1990,1995,2000,2005,2010,2015),
                                                    effect=c(summary(inpartycombined.dem.raw$model)$coefficients[16:26,1]),
                                                    se=c(summary(inpartycombined.dem.raw$model)$coefficients[16:26,2]))
appendixperiodinpartycombinedapci.dem$lci <- appendixperiodinpartycombinedapci.dem$effect - (1.96*appendixperiodinpartycombinedapci.dem$se)
appendixperiodinpartycombinedapci.dem$uci <- appendixperiodinpartycombinedapci.dem$effect + (1.96*appendixperiodinpartycombinedapci.dem$se)
appendixperiodinpartycombinedapci.dem$color <- ifelse(appendixperiodinpartycombinedapci.dem$lci<0 & appendixperiodinpartycombinedapci.dem$uci>0,"Gray50","Black")
appendixperiodinpartycombinedapci.dem$group <- 1

figurea38.a <- ggplot(appendixperiodinpartycombinedapci.dem,aes(x=period,y=effect,color=color,group=group)) + geom_point() + geom_errorbar(aes(ymin=lci,ymax=uci),size=0.5,width=0) + theme_classic() + theme(plot.title = element_text(hjust = 0.5)) + ggtitle("APCI") + geom_hline(yintercept=0,linetype="dashed",color="red") + scale_color_identity() + xlab("Year") + ylab("Deviation from Grand Mean") + ylim(-34,39)
figurea38.b <- ggplot(appendixperiodinpartycombinedie.dem,aes(x=period,y=effect,color=color,group=group)) + geom_point() + geom_errorbar(aes(ymin=lci,ymax=uci),size=0.5,width=0) + theme_classic() + theme(plot.title = element_text(hjust = 0.5)) + ggtitle("IE") + geom_hline(yintercept=0,linetype="dashed",color="red") + scale_color_identity() + xlab("Year") + ylab("Deviation from Grand Mean") + ylim(-34,39)
ggarrange(figurea38.a,figurea38.b,nrow=1)

# Figure A39
appendixcohortinpartycombinedie.dem <- data.frame(cohort=c(1870,1875,1880,1885,1890,1895,1900,1905,1910,1915,1920,1925,1930,1935,1940,1945,1950,1955,1960,1965,1970,1975,1980,1985,1990,1995),
                                                  effect=c(8.022,-4.263,-0.691,1.513,2.115,1.211,1.4446,0.959,0.131,1.029,0.258,0.504,-0.247,-1.111,-1.139,-2.412,-2.043,-1.818,-0.903,-0.450,-0.871,-0.093,-1.498,0.643,0.174,-0.465),
                                                  se=c(15.434,8.091,4.354,2.886,2.454,2.127,1.926,1.731,1.565,1.415,1.264,1.123,0.996,0.867,0.724,0.590,0.481,0.425,0.437,0.522,0.656,0.774,0.907,1.049,1.212,1.786))
appendixcohortinpartycombinedie.dem$lci <- appendixcohortinpartycombinedie.dem$effect - (1.96*appendixcohortinpartycombinedie.dem$se)
appendixcohortinpartycombinedie.dem$uci <- appendixcohortinpartycombinedie.dem$effect + (1.96*appendixcohortinpartycombinedie.dem$se)
appendixcohortinpartycombinedie.dem$color <- ifelse(appendixcohortinpartycombinedie.dem$lci<0 & appendixcohortinpartycombinedie.dem$uci>0,"Gray50","Black")
appendixcohortinpartycombinedie.dem$group <- 1

appendixcohortinpartycombinedapci.dem <- data.frame(cohort=c(1870,1875,1880,1885,1890,1895,1900,1905,1910,1915,1920,1925,1930,1935,1940,1945,1950,1955,1960,1965,1970,1975,1980,1985,1990,1995),
                                                    effect=c(as.numeric(inpartycombined.dem.raw.cohort$cohort_average$cohort_average)),
                                                    se=c(as.numeric(inpartycombined.dem.raw.cohort$cohort_average$cohort_average_se)))
appendixcohortinpartycombinedapci.dem$lci <- appendixcohortinpartycombinedapci.dem$effect - (1.96*appendixcohortinpartycombinedapci.dem$se)
appendixcohortinpartycombinedapci.dem$uci <- appendixcohortinpartycombinedapci.dem$effect + (1.96*appendixcohortinpartycombinedapci.dem$se)
appendixcohortinpartycombinedapci.dem$color <- ifelse(appendixcohortinpartycombinedapci.dem$lci<0 & appendixcohortinpartycombinedapci.dem$uci>0,"Gray50","Black")
appendixcohortinpartycombinedapci.dem$group <- 1

figurea39.a <- ggplot(appendixcohortinpartycombinedie.dem,aes(x=cohort,y=effect,color=color,group=group)) + geom_point() + geom_errorbar(aes(ymin=lci,ymax=uci),size=0.5,width=0) + theme_classic() + theme(plot.title = element_text(hjust = 0.5)) + ggtitle("APCI") + geom_hline(yintercept=0,linetype="dashed",color="red") + scale_color_identity() + xlab("Birth Cohort") + ylab("Deviation from Grand Mean") + ylim(-34,39)
figurea39.b <- ggplot(appendixcohortinpartycombinedapci.dem,aes(x=cohort,y=effect,color=color,group=group)) + geom_point() + geom_errorbar(aes(ymin=lci,ymax=uci),size=0.5,width=0) + theme_classic() + theme(plot.title = element_text(hjust = 0.5)) + ggtitle("IE") + geom_hline(yintercept=0,linetype="dashed",color="red") + scale_color_identity() + xlab("Birth Cohort") + ylab("Deviation from Grand Mean") + ylim(-34,39)
ggarrange(figurea39.a,figurea39.b,nrow=1)

# Figure A40
appendixcohortinpartycombinedapcislopes.dem <- data.frame(cohort=c(1870,1875,1880,1885,1890,1895,1900,1905,1910,1915,1920,1925,1930,1935,1940,1945,1950,1955,1960,1965,1970,1975,1980,1985,1990,1995),
                                                          effect=c(as.numeric(inpartycombined.dem.raw.cohort$cohort_slope$cohort_slope)),
                                                          se=c(as.numeric(inpartycombined.dem.raw.cohort$cohort_slope$cohort_slope_se)))
appendixcohortinpartycombinedapcislopes.dem$lci <- appendixcohortinpartycombinedapcislopes.dem$effect - (1.96*appendixcohortinpartycombinedapcislopes.dem$se)
appendixcohortinpartycombinedapcislopes.dem$uci <- appendixcohortinpartycombinedapcislopes.dem$effect + (1.96*appendixcohortinpartycombinedapcislopes.dem$se)
appendixcohortinpartycombinedapcislopes.dem$color <- ifelse(appendixcohortinpartycombinedapcislopes.dem$lci<0 & appendixcohortinpartycombinedapcislopes.dem$uci>0,"Gray50","Black")
appendixcohortinpartycombinedapcislopes.dem$group <- 1

ggplot(appendixcohortinpartycombinedapcislopes.dem,aes(x=cohort,y=effect,color=color,group=group)) + geom_point() + geom_errorbar(aes(ymin=lci,ymax=uci),size=0.5,width=0) + theme_classic() + theme(plot.title = element_text(hjust = 0.5)) + geom_hline(yintercept=0,linetype="dashed",color="red") + scale_color_identity() + xlab("Birth Cohort") + ylab("Deviation from Grand Mean") + ylim(-34,39)

# Figure A41
appendixageinpartycombinedie.rep <- data.frame(age=c(20,25,30,35,40,45,50,55,60,65,70,75,80,85,90),
                                               effect=c(-2.564,-2.637,-2.886,-2.325,-0.866,-1.471,-1.271,-0.185,0.912,1.520,1.629,1.501,2.276,1.001,5.365),
                                               se=c(0.82,0.719,0.639,0.577,0.542,0.525,0.508,0.522,0.558,0.619,0.713,0.827,1.015,1.417,2.235))
appendixageinpartycombinedie.rep$lci <- appendixageinpartycombinedie.rep$effect - (1.96*appendixageinpartycombinedie.rep$se)
appendixageinpartycombinedie.rep$uci <- appendixageinpartycombinedie.rep$effect + (1.96*appendixageinpartycombinedie.rep$se)
appendixageinpartycombinedie.rep$color <- ifelse(appendixageinpartycombinedie.rep$lci<0 & appendixageinpartycombinedie.rep$uci>0,"Gray50","Black")
appendixageinpartycombinedie.rep$group <- 1

appendixageinpartycombinedapci.rep <- data.frame(age=c(25,30,35,40,45,50,55,60,65,70,75,80,85),
                                                 effect=c(summary(inpartycombined.rep.raw$model)$coefficients[2:14,1]),
                                                 se=c(summary(inpartycombined.rep.raw$model)$coefficients[2:14,2]))
appendixageinpartycombinedapci.rep$lci <- appendixageinpartycombinedapci.rep$effect - (1.96*appendixageinpartycombinedapci.rep$se)
appendixageinpartycombinedapci.rep$uci <- appendixageinpartycombinedapci.rep$effect + (1.96*appendixageinpartycombinedapci.rep$se)
appendixageinpartycombinedapci.rep$color <- ifelse(appendixageinpartycombinedapci.rep$lci<0 & appendixageinpartycombinedapci.rep$uci>0,"Gray50","Black")
appendixageinpartycombinedapci.rep$group <- 1

figurea41.a <- ggplot(appendixageinpartycombinedapci.rep,aes(x=age,y=effect,color=color,group=group)) + geom_point() + geom_errorbar(aes(ymin=lci,ymax=uci),size=0.5,width=0) + theme_classic() + theme(plot.title = element_text(hjust = 0.5)) + ggtitle("APCI") + geom_hline(yintercept=0,linetype="dashed",color="red") + scale_color_identity() + xlab("Age") + ylab("Deviation from Grand Mean") + ylim(-34,39)
figurea41.b <- ggplot(appendixageinpartycombinedie.rep,aes(x=age,y=effect,color=color,group=group)) + geom_point() + geom_errorbar(aes(ymin=lci,ymax=uci),size=0.5,width=0) + theme_classic() + theme(plot.title = element_text(hjust = 0.5)) + ggtitle("IE") + geom_hline(yintercept=0,linetype="dashed",color="red") + scale_color_identity() + xlab("Age") + ylab("Deviation from Grand Mean") + ylim(-34,39)
ggarrange(figurea41.a,figurea41.b,nrow=1)

## Figure A42
appendixperiodinpartycombinedie.rep <- data.frame(period=c(1960,1965,1970,1975,1980,1985,1990,1995,2000,2005,2010,2015),
                                                  effect=c(6.689,4.84,-0.306,-2.244,2.773,2.552,-1.453,-0.763,3.298,-3.261,-3.967,-8.158),
                                                  se=c(0.886,0.647,0.608,0.497,0.509,0.389,0.434,0.441,0.789,0.739,0.603,0.710))
appendixperiodinpartycombinedie.rep$lci <- appendixperiodinpartycombinedie.rep$effect - (1.96*appendixperiodinpartycombinedie.rep$se)
appendixperiodinpartycombinedie.rep$uci <- appendixperiodinpartycombinedie.rep$effect + (1.96*appendixperiodinpartycombinedie.rep$se)
appendixperiodinpartycombinedie.rep$color <- ifelse(appendixperiodinpartycombinedie.rep$lci<0 & appendixperiodinpartycombinedie.rep$uci>0,"Gray50","Black")
appendixperiodinpartycombinedie.rep$group <- 1

appendixperiodinpartycombinedapci.rep <- data.frame(period=c(1965,1970,1975,1980,1985,1990,1995,2000,2005,2010,2015),
                                                    effect=c(summary(inpartycombined.rep.raw$model)$coefficients[15:25,1]),
                                                    se=c(summary(inpartycombined.rep.raw$model)$coefficients[15:25,2]))
appendixperiodinpartycombinedapci.rep$lci <- appendixperiodinpartycombinedapci.rep$effect - (1.96*appendixperiodinpartycombinedapci.rep$se)
appendixperiodinpartycombinedapci.rep$uci <- appendixperiodinpartycombinedapci.rep$effect + (1.96*appendixperiodinpartycombinedapci.rep$se)
appendixperiodinpartycombinedapci.rep$color <- ifelse(appendixperiodinpartycombinedapci.rep$lci<0 & appendixperiodinpartycombinedapci.rep$uci>0,"Gray50","Black")
appendixperiodinpartycombinedapci.rep$group <- 1

figurea42.a <- ggplot(appendixperiodinpartycombinedapci.rep,aes(x=period,y=effect,color=color,group=group)) + geom_point() + geom_errorbar(aes(ymin=lci,ymax=uci),size=0.5,width=0) + theme_classic() + theme(plot.title = element_text(hjust = 0.5)) + ggtitle("APCI") + geom_hline(yintercept=0,linetype="dashed",color="red") + scale_color_identity() + xlab("Year") + ylab("Deviation from Grand Mean") + ylim(-34,39)
figurea42.b <- ggplot(appendixperiodinpartycombinedie.rep,aes(x=period,y=effect,color=color,group=group)) + geom_point() + geom_errorbar(aes(ymin=lci,ymax=uci),size=0.5,width=0) + theme_classic() + theme(plot.title = element_text(hjust = 0.5)) + ggtitle("IE") + geom_hline(yintercept=0,linetype="dashed",color="red") + scale_color_identity() + xlab("Year") + ylab("Deviation from Grand Mean") + ylim(-34,39)
ggarrange(figurea42.a,figurea42.b,nrow=1)

# Figure A43
appendixcohortinpartycombinedie.rep <- data.frame(cohort=c(1870,1875,1880,1885,1890,1895,1900,1905,1910,1915,1920,1925,1930,1935,1940,1945,1950,1955,1960,1965,1970,1975,1980,1985,1990,1995),
                                                  effect=c(5.548,-13.463,7.914,3.153,4.001,3.697,0.639,1.845,0.964,0.633,-1.151,-1.335,-0.196,-0.510,-1.975,-1.622,-1.946,-1.419,-1.458,-2.048,-0.459,-2.194,-0.625,-1.867,1.555,2.319),
                                                  se=c(9.092,7.132,4.224,2.504,1.999,1.709,1.569,1.405,1.260,1.146,1.052,0.930,0.863,0.759,0.663,0.58,0.528,0.483,0.493,0.573,0.7,0.854,1.025,1.056,1.22,1.871))
appendixcohortinpartycombinedie.rep$lci <- appendixcohortinpartycombinedie.rep$effect - (1.96*appendixcohortinpartycombinedie.rep$se)
appendixcohortinpartycombinedie.rep$uci <- appendixcohortinpartycombinedie.rep$effect + (1.96*appendixcohortinpartycombinedie.rep$se)
appendixcohortinpartycombinedie.rep$color <- ifelse(appendixcohortinpartycombinedie.rep$lci<0 & appendixcohortinpartycombinedie.rep$uci>0,"Gray50","Black")
appendixcohortinpartycombinedie.rep$group <- 1

appendixcohortinpartycombinedapci.rep <- data.frame(cohort=c(1875,1880,1885,1890,1895,1900,1905,1910,1915,1920,1925,1930,1935,1940,1945,1950,1955,1960,1965,1970,1975,1980,1985,1990,1995),
                                                    effect=c(as.numeric(inpartycombined.rep.raw.cohort$cohort_average$cohort_average)),
                                                    se=c(as.numeric(inpartycombined.rep.raw.cohort$cohort_average$cohort_average_se)))
appendixcohortinpartycombinedapci.rep$lci <- appendixcohortinpartycombinedapci.rep$effect - (1.96*appendixcohortinpartycombinedapci.rep$se)
appendixcohortinpartycombinedapci.rep$uci <- appendixcohortinpartycombinedapci.rep$effect + (1.96*appendixcohortinpartycombinedapci.rep$se)
appendixcohortinpartycombinedapci.rep$color <- ifelse(appendixcohortinpartycombinedapci.rep$lci<0 & appendixcohortinpartycombinedapci.rep$uci>0,"Gray50","Black")
appendixcohortinpartycombinedapci.rep$group <- 1

figurea43.a <- ggplot(appendixcohortinpartycombinedapci.rep,aes(x=cohort,y=effect,color=color,group=group)) + geom_point() + geom_errorbar(aes(ymin=lci,ymax=uci),size=0.5,width=0) + theme_classic() + theme(plot.title = element_text(hjust = 0.5)) + ggtitle("APCI") + geom_hline(yintercept=0,linetype="dashed",color="red") + scale_color_identity() + xlab("Birth Cohort") + ylab("Deviation from Grand Mean") + ylim(-34,39)
figurea43.b <- ggplot(appendixcohortinpartycombinedie.rep,aes(x=cohort,y=effect,color=color,group=group)) + geom_point() + geom_errorbar(aes(ymin=lci,ymax=uci),size=0.5,width=0) + theme_classic() + theme(plot.title = element_text(hjust = 0.5)) + ggtitle("IE") + geom_hline(yintercept=0,linetype="dashed",color="red") + scale_color_identity() + xlab("Birth Cohort") + ylab("Deviation from Grand Mean") + ylim(-34,39)
ggarrange(figurea43.a,figurea43.b,nrow=1)

# Figure A44
appendixcohortinpartycombinedapcislopes.rep <- data.frame(cohort=c(1875,1880,1885,1890,1895,1900,1905,1910,1915,1920,1925,1930,1935,1940,1945,1950,1955,1960,1965,1970,1975,1980,1985,1990,1995),
                                                          effect=c(as.numeric(inpartycombined.rep.raw.cohort$cohort_slope$cohort_slope)),
                                                          se=c(as.numeric(inpartycombined.rep.raw.cohort$cohort_slope$cohort_slope_se)))
appendixcohortinpartycombinedapcislopes.rep$lci <- appendixcohortinpartycombinedapcislopes.rep$effect - (1.96*appendixcohortinpartycombinedapcislopes.rep$se)
appendixcohortinpartycombinedapcislopes.rep$uci <- appendixcohortinpartycombinedapcislopes.rep$effect + (1.96*appendixcohortinpartycombinedapcislopes.rep$se)
appendixcohortinpartycombinedapcislopes.rep$color <- ifelse(appendixcohortinpartycombinedapcislopes.rep$lci<0 & appendixcohortinpartycombinedapcislopes.rep$uci>0,"Gray50","Black")
appendixcohortinpartycombinedapcislopes.rep$group <- 1

ggplot(appendixcohortinpartycombinedie.rep,aes(x=cohort,y=effect,color=color,group=group)) + geom_point() + geom_errorbar(aes(ymin=lci,ymax=uci),size=0.5,width=0) + theme_classic() + theme(plot.title = element_text(hjust = 0.5)) + geom_hline(yintercept=0,linetype="dashed",color="red") + scale_color_identity() + xlab("Birth Cohort") + ylab("Deviation from Grand Mean") + ylim(-34,39)

# Table A27
## Model 1
outpartynew.all.raw <- APCI::temp_model(data=anes.cumulative.full,outcome='outparty.new',age='agegroup',period='yeargroup',cohort='cohort',covariate=c(),family='gaussian')
outpartynew.all.raw.cohort <- APCI::cohortdeviation(outpartynew.all.raw$A,outpartynew.all.raw$P,outpartynew.all.raw$C,model=outpartynew.all.raw$model,covariate=c())
## Model 2
outpartynewpid.all.raw <- APCI::temp_model(data=anes.cumulative.full,outcome='outparty.new',age='agegroup',period='yeargroup',cohort='cohort',covariate=c('strong'),family='gaussian')
outpartynewpid.all.raw.cohort <- APCI::cohortdeviation(outpartynewpid.all.raw$A,outpartynewpid.all.raw$P,outpartynewpid.all.raw$C,model=outpartynewpid.all.raw$model,covariate=c('strong'))
## Model 3
outpartynewsort.all.raw <- APCI::temp_model(data=anes.cumulative.full,outcome='outparty.new',age='agegroup',period='yeargroup',cohort='cohort',covariate=c('ideosort','socialsort'),family='gaussian')
outpartynewsort.all.raw.cohort <- APCI::cohortdeviation(outpartynewsort.all.raw$A,outpartynewsort.all.raw$P,outpartynewsort.all.raw$C,model=outpartynewsort.all.raw$model,covariate=c('strong'))
## Model 4
outpartynew.all.controls <- APCI::temp_model(data=anes.cumulative.full,outcome='outparty.new',age='agegroup',period='yeargroup',cohort='cohort',covariate=c('gender','black','education','ideo','attendance','south'),family='gaussian')
outpartynew.all.controls.cohort <- APCI::cohortdeviation(outpartynew.all.controls$A,outpartynew.all.controls$P,outpartynew.all.controls$C,model=outpartynew.all.controls$model,covariate=c('gender','black','education','ideo','attendance','south'))
## Model 5
outpartynewpid.all.controls <- APCI::temp_model(data=anes.cumulative.full,outcome='outparty.new',age='agegroup',period='yeargroup',cohort='cohort',covariate=c('strong','gender','black','education','ideo','attendance','south'),family='gaussian')
outpartynewpid.all.controls.cohort <- APCI::cohortdeviation(outpartynewpid.all.controls$A,outpartynewpid.all.controls$P,outpartynewpid.all.controls$C,model=outpartynewpid.all.controls$model,covariate=c('strong','gender','black','education','ideo','attendance','south'))
## Model 6
outpartynewsort.all.controls <- APCI::temp_model(data=anes.cumulative.full,outcome='outparty.new',age='agegroup',period='yeargroup',cohort='cohort',covariate=c('ideosort','socialsort','gender','black','education','ideo','attendance','south'),family='gaussian')
outpartynewsort.all.controls.cohort <- APCI::cohortdeviation(outpartynewsort.all.controls$A,outpartynewsort.all.controls$P,outpartynewsort.all.controls$C,model=outpartynewsort.all.controls$model,covariate=c('ideosort','socialsort','gender','black','education','ideo','attendance','south'))

# Table A29
## Model 1
outpartynew.dem.raw <- APCI::temp_model(data=anes.cumulative.dem,outcome='outparty.new',age='agegroup',period='yeargroup',cohort='cohort',covariate=c(),family='gaussian')
outpartynew.dem.raw.cohort <- APCI::cohortdeviation(outpartynew.dem.raw$A,outpartynew.dem.raw$P,outpartynew.dem.raw$C,model=outpartynew.dem.raw$model,covariate=c())
## Model 2
outpartynewpid.dem.raw <- APCI::temp_model(data=anes.cumulative.dem,outcome='outparty.new',age='agegroup',period='yeargroup',cohort='cohort',covariate=c('strong'),family='gaussian')
outpartynewpid.dem.raw.cohort <- APCI::cohortdeviation(outpartynewpid.dem.raw$A,outpartynewpid.dem.raw$P,outpartynewpid.dem.raw$C,model=outpartynewpid.dem.raw$model,covariate=c('strong'))
## Model 3
outpartynewsort.dem.raw <- APCI::temp_model(data=anes.cumulative.dem,outcome='outparty.new',age='agegroup',period='yeargroup',cohort='cohort',covariate=c('ideosort','socialsort'),family='gaussian')
outpartynewsort.dem.raw.cohort <- APCI::cohortdeviation(outpartynewsort.dem.raw$A,outpartynewsort.dem.raw$P,outpartynewsort.dem.raw$C,model=outpartynewsort.dem.raw$model,covariate=c('strong'))
## Model 4
outpartynew.dem.controls <- APCI::temp_model(data=anes.cumulative.dem,outcome='outparty.new',age='agegroup',period='yeargroup',cohort='cohort',covariate=c('gender','black','education','ideo','attendance','south'),family='gaussian')
outpartynew.dem.controls.cohort <- APCI::cohortdeviation(outpartynew.dem.controls$A,outpartynew.dem.controls$P,outpartynew.dem.controls$C,model=outpartynew.dem.controls$model,covariate=c('gender','black','education','ideo','attendance','south'))
## Model 5
outpartynewpid.dem.controls <- APCI::temp_model(data=anes.cumulative.dem,outcome='outparty.new',age='agegroup',period='yeargroup',cohort='cohort',covariate=c('strong','gender','black','education','ideo','attendance','south'),family='gaussian')
outpartynewpid.dem.controls.cohort <- APCI::cohortdeviation(outpartynewpid.dem.controls$A,outpartynewpid.dem.controls$P,outpartynewpid.dem.controls$C,model=outpartynewpid.dem.controls$model,covariate=c('strong','gender','black','education','ideo','attendance','south'))
## Model 6
outpartynewsort.dem.controls <- APCI::temp_model(data=anes.cumulative.dem,outcome='outparty.new',age='agegroup',period='yeargroup',cohort='cohort',covariate=c('ideosort','socialsort','gender','black','education','ideo','attendance','south'),family='gaussian')
outpartynewsort.dem.controls.cohort <- APCI::cohortdeviation(outpartynewsort.dem.controls$A,outpartynewsort.dem.controls$P,outpartynewsort.dem.controls$C,model=outpartynewsort.dem.controls$model,covariate=c('ideosort','socialsort','gender','black','education','ideo','attendance','south'))

# Table A31
## Model 1
outpartynew.rep.raw <- APCI::temp_model(data=repnot90,outcome='outparty.new',age='agegroup',period='yeargroup',cohort='cohort',covariate=c(),family='gaussian')
outpartynew.rep.raw.cohort <- APCI::cohortdeviation(outpartynew.rep.raw$A,outpartynew.rep.raw$P,outpartynew.rep.raw$C,model=outpartynew.rep.raw$model,covariate=c())
## Model 2
outpartynewpid.rep.raw <- APCI::temp_model(data=repnot90,outcome='outparty.new',age='agegroup',period='yeargroup',cohort='cohort',covariate=c('strong'),family='gaussian')
outpartynewpid.rep.raw.cohort <- APCI::cohortdeviation(outpartynewpid.rep.raw$A,outpartynewpid.rep.raw$P,outpartynewpid.rep.raw$C,model=outpartynewpid.rep.raw$model,covariate=c('strong'))
## Model 3
outpartynewsort.rep.raw <- APCI::temp_model(data=repnot90,outcome='outparty.new',age='agegroup',period='yeargroup',cohort='cohort',covariate=c('ideosort','socialsort'),family='gaussian')
outpartynewsort.rep.raw.cohort <- APCI::cohortdeviation(outpartynewsort.rep.raw$A,outpartynewsort.rep.raw$P,outpartynewsort.rep.raw$C,model=outpartynewsort.rep.raw$model,covariate=c('strong'))
## Model 4
outpartynew.rep.controls <- APCI::temp_model(data=repnot90,outcome='outparty.new',age='agegroup',period='yeargroup',cohort='cohort',covariate=c('gender','black','education','ideo','attendance','south'),family='gaussian')
outpartynew.rep.controls.cohort <- APCI::cohortdeviation(outpartynew.rep.controls$A,outpartynew.rep.controls$P,outpartynew.rep.controls$C,model=outpartynew.rep.controls$model,covariate=c('gender','black','education','ideo','attendance','south'))
## Model 5
outpartynewpid.rep.controls <- APCI::temp_model(data=repnot90,outcome='outparty.new',age='agegroup',period='yeargroup',cohort='cohort',covariate=c('strong','gender','black','education','ideo','attendance','south'),family='gaussian')
outpartynewpid.rep.controls.cohort <- APCI::cohortdeviation(outpartynewpid.rep.controls$A,outpartynewpid.rep.controls$P,outpartynewpid.rep.controls$C,model=outpartynewpid.rep.controls$model,covariate=c('strong','gender','black','education','ideo','attendance','south'))
## Model 6
outpartynewsort.rep.controls <- APCI::temp_model(data=repnot90,outcome='outparty.new',age='agegroup',period='yeargroup',cohort='cohort',covariate=c('ideosort','socialsort','gender','black','education','ideo','attendance','south'),family='gaussian')
outpartynewsort.rep.controls.cohort <- APCI::cohortdeviation(outpartynewsort.rep.controls$A,outpartynewsort.rep.controls$P,outpartynewsort.rep.controls$C,model=outpartynewsort.rep.controls$model,covariate=c('ideosort','socialsort','gender','black','education','ideo','attendance','south'))

# Figure A45
appendixageoutpartynewie.all <- data.frame(age=c(20,25,30,35,40,45,50,55,60,65,70,75,80,85,90),
                                           effect=c(2.870,1.496,2.307,1.164,0.517,0.001,0.871,-0.477,-0.402,-1.453,-1.593,-0.522,-0.437,-1.984,-2.358),
                                           se=c(0.694,0.614,0.546,0.496,0.482,0.474,0.465,0.478,0.519,0.565,0.658,0.767,0.92,1.25,2.215))
appendixageoutpartynewie.all$lci <- appendixageoutpartynewie.all$effect - (1.96*appendixageoutpartynewie.all$se)
appendixageoutpartynewie.all$uci <- appendixageoutpartynewie.all$effect + (1.96*appendixageoutpartynewie.all$se)
appendixageoutpartynewie.all$color <- ifelse(appendixageoutpartynewie.all$lci<0 & appendixageoutpartynewie.all$uci>0,"Gray50","Black")
appendixageoutpartynewie.all$group <- 1

appendixageoutpartynewapci.all <- data.frame(age=c(25,30,35,40,45,50,55,60,65,70,75,80,85,90),
                                             effect=c(summary(outpartynew.all.raw$model)$coefficients[2:15,1]),
                                             se=c(summary(outpartynew.all.raw$model)$coefficients[2:15,1]))
appendixageoutpartynewapci.all$lci <- appendixageoutpartynewapci.all$effect - (1.96*appendixageoutpartynewapci.all$se)
appendixageoutpartynewapci.all$uci <- appendixageoutpartynewapci.all$effect + (1.96*appendixageoutpartynewapci.all$se)
appendixageoutpartynewapci.all$color <- ifelse(appendixageoutpartynewapci.all$lci<0 & appendixageoutpartynewapci.all$uci>0,"Gray50","Black")
appendixageoutpartynewapci.all$group <- 1

figurea45.a <- ggplot(appendixageoutpartynewapci.all,aes(x=age,y=effect,color=color,group=group)) + geom_point() + geom_errorbar(aes(ymin=lci,ymax=uci),size=0.5,width=0) + theme_classic() + theme(plot.title = element_text(hjust = 0.5)) + ggtitle("APCI") + geom_hline(yintercept=0,linetype="dashed",color="red") + scale_color_identity() + xlab("Age") + ylab("Deviation from Grand Mean") + ylim(-63,28)
figurea45.b <- ggplot(appendixageoutpartynewie.all,aes(x=age,y=effect,color=color,group=group)) + geom_point() + geom_errorbar(aes(ymin=lci,ymax=uci),size=0.5,width=0) + theme_classic() + theme(plot.title = element_text(hjust = 0.5)) + ggtitle("IE") + geom_hline(yintercept=0,linetype="dashed",color="red") + scale_color_identity() + xlab("Age") + ylab("Deviation from Grand Mean") + ylim(-63,28)
ggarrange(figurea45.a,figurea45.b,nrow=1)

# Figure A46
appendixperiodoutpartynewie.all <- data.frame(period=c(1975,1980,1985,1990,1995,2000,2005,2010,2015),
                                              effect=c(7.763,5.896,6.991,2.767,2.15,0.323,-2.983,-10.695,-12.212),
                                              se=c(0.483,0.45,0.35,0.357,0.326,0.601,0.479,0.401,0.503))
appendixperiodoutpartynewie.all$lci <- appendixperiodoutpartynewie.all$effect - (1.96*appendixperiodoutpartynewie.all$se)
appendixperiodoutpartynewie.all$uci <- appendixperiodoutpartynewie.all$effect + (1.96*appendixperiodoutpartynewie.all$se)
appendixperiodoutpartynewie.all$color <- ifelse(appendixperiodoutpartynewie.all$lci<0 & appendixperiodoutpartynewie.all$uci>0,"Gray50","Black")
appendixperiodoutpartynewie.all$group <- 1

appendixperiodoutpartynewapci.all <- data.frame(period=c(1980,1985,1990,1995,2000,2005,2010,2015),
                                                effect=c(summary(outpartynew.all.raw$model)$coefficients[16:23,1]),
                                                se=c(summary(outpartynew.all.raw$model)$coefficients[16:23,2]))
appendixperiodoutpartynewapci.all$lci <- appendixperiodoutpartynewapci.all$effect - (1.96*appendixperiodoutpartynewapci.all$se)
appendixperiodoutpartynewapci.all$uci <- appendixperiodoutpartynewapci.all$effect + (1.96*appendixperiodoutpartynewapci.all$se)
appendixperiodoutpartynewapci.all$color <- ifelse(appendixperiodoutpartynewapci.all$lci<0 & appendixperiodoutpartynewapci.all$uci>0,"Gray50","Black")
appendixperiodoutpartynewapci.all$group <- 1

figurea46.a <- ggplot(appendixperiodoutpartynewapci.all,aes(x=period,y=effect,color=color,group=group)) + geom_point() + geom_errorbar(aes(ymin=lci,ymax=uci),size=0.5,width=0) + theme_classic() + theme(plot.title = element_text(hjust = 0.5)) + ggtitle("APCI") + geom_hline(yintercept=0,linetype="dashed",color="red") + scale_color_identity() + xlab("Year") + ylab("Deviation from Grand Mean") + ylim(-63,28)
figurea46.b <- ggplot(appendixperiodoutpartynewie.all,aes(x=period,y=effect,color=color,group=group)) + geom_point() + geom_errorbar(aes(ymin=lci,ymax=uci),size=0.5,width=0) + theme_classic() + theme(plot.title = element_text(hjust = 0.5)) + ggtitle("IE") + geom_hline(yintercept=0,linetype="dashed",color="red") + scale_color_identity() + xlab("Year") + ylab("Deviation from Grand Mean") + ylim(-63,28)
ggarrange(figurea46.a,figurea46.b,nrow=1)

# Figure A47
appendixcohortoutpartynewie.all <- data.frame(cohort=c(1885,1890,1895,1900,1905,1910,1915,1920,1925,1930,1935,1940,1945,1950,1955,1960,1965,1970,1975,1980,1985,1990,1995),
                                              effect=c(-19.498,15.444,-1.515,-0.532,1.481,2.839,2.082,2.784,1.847,0.648,1.835,0.107,-1.131,-1.198,-2.300,-1.267,-0.484,0.167,-1.583,-1.55,-0.216,-0.926,2.964),
                                              se=c(6.788,4.171,2.322,1.837,1.413,1.227,1.113,1.006,0.902,0.844,0.757,0.657,0.561,0.483,0.42,0.409,0.455,0.542,0.625,0.721,0.783,0.9,1.457))
appendixcohortoutpartynewie.all$lci <- appendixcohortoutpartynewie.all$effect - (1.96*appendixcohortoutpartynewie.all$se)
appendixcohortoutpartynewie.all$uci <- appendixcohortoutpartynewie.all$effect + (1.96*appendixcohortoutpartynewie.all$se)
appendixcohortoutpartynewie.all$color <- ifelse(appendixcohortoutpartynewie.all$lci<0 & appendixcohortoutpartynewie.all$uci>0,"Gray50","Black")
appendixcohortoutpartynewie.all$group <- 1

appendixcohortoutpartynewapci.all <- data.frame(cohort=c(1885,1890,1895,1900,1905,1910,1915,1920,1925,1930,1935,1940,1945,1950,1955,1960,1965,1970,1975,1980,1985,1990,1995),
                                                effect=c(as.numeric(outpartynew.all.raw.cohort$cohort_average$cohort_average)),
                                                se=c(as.numeric(outpartynew.all.raw.cohort$cohort_average$cohort_average_se)))
appendixcohortoutpartynewapci.all$lci <- appendixcohortoutpartynewapci.all$effect - (1.96*appendixcohortoutpartynewapci.all$se)
appendixcohortoutpartynewapci.all$uci <- appendixcohortoutpartynewapci.all$effect + (1.96*appendixcohortoutpartynewapci.all$se)
appendixcohortoutpartynewapci.all$color <- ifelse(appendixcohortoutpartynewapci.all$lci<0 & appendixcohortoutpartynewapci.all$uci>0,"Gray50","Black")
appendixcohortoutpartynewapci.all$group <- 1

figurea47.a <- ggplot(appendixcohortoutpartynewapci.all,aes(x=cohort,y=effect,color=color,group=group)) + geom_point() + geom_errorbar(aes(ymin=lci,ymax=uci),size=0.5,width=0) + theme_classic() + theme(plot.title = element_text(hjust = 0.5)) + ggtitle("APCI") + geom_hline(yintercept=0,linetype="dashed",color="red") + scale_color_identity() + xlab("Birth Cohort") + ylab("Deviation from Grand Mean") + ylim(-63,28)
figurea47.b <- ggplot(appendixcohortoutpartynewie.all,aes(x=cohort,y=effect,color=color,group=group)) + geom_point() + geom_errorbar(aes(ymin=lci,ymax=uci),size=0.5,width=0) + theme_classic() + theme(plot.title = element_text(hjust = 0.5)) + ggtitle("IE") + geom_hline(yintercept=0,linetype="dashed",color="red") + scale_color_identity() + xlab("Birth Cohort") + ylab("Deviation from Grand Mean") + ylim(-63,28)
ggarrange(figurea47.a,figurea47.b,nrow=1)

# Figure A48
appendixcohortoutpartynewapcislopes.all <- data.frame(cohort=c(1885,1890,1895,1900,1905,1910,1915,1920,1925,1930,1935,1940,1945,1950,1955,1960,1965,1970,1975,1980,1985,1990,1995),
                                                      effect=c(as.numeric(outpartynew.all.raw.cohort$cohort_slope$cohort_slope)),
                                                      se=c(as.numeric(outpartynew.all.raw.cohort$cohort_slope$cohort_slope_se)))
appendixcohortoutpartynewapcislopes.all$lci <- appendixcohortoutpartynewapcislopes.all$effect - (1.96*appendixcohortoutpartynewapcislopes.all$se)
appendixcohortoutpartynewapcislopes.all$uci <- appendixcohortoutpartynewapcislopes.all$effect + (1.96*appendixcohortoutpartynewapcislopes.all$se)
appendixcohortoutpartynewapcislopes.all$color <- ifelse(appendixcohortoutpartynewapcislopes.all$lci<0 & appendixcohortoutpartynewapcislopes.all$uci>0,"Gray50","Black")
appendixcohortoutpartynewapcislopes.all$group <- 1

ggplot(appendixcohortoutpartynewapcislopes.all,aes(x=cohort,y=effect,color=color,group=group)) + geom_point() + geom_errorbar(aes(ymin=lci,ymax=uci),size=0.5,width=0) + theme_classic() + theme(plot.title = element_text(hjust = 0.5)) + geom_hline(yintercept=0,linetype="dashed",color="red") + scale_color_identity() + xlab("Birth Cohort") + ylab("Deviation from Grand Mean") + ylim(-63,28)

# Figure A49
appendixageoutpartynewie.dem <- data.frame(age=c(20,25,30,35,40,45,50,55,60,65,70,75,80,85,90),
                                           effect=c(1.878,0.967,1.78,0.911,0.007,-0.908,0.865,-0.383,0.036,-0.994,-0.452,1.760,0.945,-1.444,-5.166),
                                           se=c(0.957,0.845,0.753,0.683,0.663,0.653,0.644,0.665,0.72,0.782,0.917,1.071,1.283,1.752,3.34))
appendixageoutpartynewie.dem$lci <- appendixageoutpartynewie.dem$effect - (1.96*appendixageoutpartynewie.dem$se)
appendixageoutpartynewie.dem$uci <- appendixageoutpartynewie.dem$effect + (1.96*appendixageoutpartynewie.dem$se)
appendixageoutpartynewie.dem$color <- ifelse(appendixageoutpartynewie.dem$lci<0 & appendixageoutpartynewie.dem$uci>0,"Gray50","Black")
appendixageoutpartynewie.dem$group <- 1

appendixageoutpartynewapci.dem <- data.frame(age=c(25,30,35,40,45,50,55,60,65,70,75,80,85,90),
                                             effect=c(summary(outpartynew.dem.raw$model)$coefficients[2:15,1]),
                                             se=c(summary(outpartynew.dem.raw$model)$coefficients[2:15,2]))
appendixageoutpartynewapci.dem$lci <- appendixageoutpartynewapci.dem$effect - (1.96*appendixageoutpartynewapci.dem$se)
appendixageoutpartynewapci.dem$uci <- appendixageoutpartynewapci.dem$effect + (1.96*appendixageoutpartynewapci.dem$se)
appendixageoutpartynewapci.dem$color <- ifelse(appendixageoutpartynewapci.dem$lci<0 & appendixageoutpartynewapci.dem$uci>0,"Gray50","Black")
appendixageoutpartynewapci.dem$group <- 1

figurea49.a <- ggplot(appendixageoutpartynewapci.dem,aes(x=age,y=effect,color=color,group=group)) + geom_point() + geom_errorbar(aes(ymin=lci,ymax=uci),size=0.5,width=0) + theme_classic() + theme(plot.title = element_text(hjust = 0.5)) + ggtitle("APCI") + geom_hline(yintercept=0,linetype="dashed",color="red") + scale_color_identity() + xlab("Age") + ylab("Deviation from Grand Mean") + ylim(-63,28)
figurea49.b <- ggplot(appendixageoutpartynewie.dem,aes(x=age,y=effect,color=color,group=group)) + geom_point() + geom_errorbar(aes(ymin=lci,ymax=uci),size=0.5,width=0) + theme_classic() + theme(plot.title = element_text(hjust = 0.5)) + ggtitle("IE") + geom_hline(yintercept=0,linetype="dashed",color="red") + scale_color_identity() + xlab("Age") + ylab("Deviation from Grand Mean") + ylim(-63,28)
ggarrange(figurea49.a,figurea49.b,nrow=1)

# Figure A50
appendixperiodoutpartynewie.dem <- data.frame(period=c(1975,1980,1985,1990,1995,2000,2005,2010,2015),
                                              effect=c(9.505,5.626,7.004,4.098,2.191,-1.155,-4.68,-11.504,-11.086),
                                              se=c(0.658,0.619,0.482,0.497,0.447,0.841,0.618,0.551,0.715))
appendixperiodoutpartynewie.dem$lci <- appendixperiodoutpartynewie.dem$effect - (1.96*appendixperiodoutpartynewie.dem$se)
appendixperiodoutpartynewie.dem$uci <- appendixperiodoutpartynewie.dem$effect + (1.96*appendixperiodoutpartynewie.dem$se)
appendixperiodoutpartynewie.dem$color <- ifelse(appendixperiodoutpartynewie.dem$lci<0 & appendixperiodoutpartynewie.dem$uci>0,"Gray50","Black")
appendixperiodoutpartynewie.dem$group <- 1

appendixperiodoutpartynewapci.dem <- data.frame(period=c(1980,1985,1990,1995,2000,2005,2010,2015),
                                                effect=c(summary(outpartynew.dem.raw$model)$coefficients[16:23,1]),
                                                se=c(summary(outpartynew.dem.raw$model)$coefficients[16:23,2]))
appendixperiodoutpartynewapci.dem$lci <- appendixperiodoutpartynewapci.dem$effect - (1.96*appendixperiodoutpartynewapci.dem$se)
appendixperiodoutpartynewapci.dem$uci <- appendixperiodoutpartynewapci.dem$effect + (1.96*appendixperiodoutpartynewapci.dem$se)
appendixperiodoutpartynewapci.dem$color <- ifelse(appendixperiodoutpartynewapci.dem$lci<0 & appendixperiodoutpartynewapci.dem$uci>0,"Gray50","Black")
appendixperiodoutpartynewapci.dem$group <- 1

figurea50.a <- ggplot(appendixperiodoutpartynewapci.dem,aes(x=period,y=effect,color=color,group=group)) + geom_point() + geom_errorbar(aes(ymin=lci,ymax=uci),size=0.5,width=0) + theme_classic() + theme(plot.title = element_text(hjust = 0.5)) + ggtitle("APCI") + geom_hline(yintercept=0,linetype="dashed",color="red") + scale_color_identity() + xlab("Year") + ylab("Deviation from Grand Mean") + ylim(-63,28)
figurea50.b <- ggplot(appendixperiodoutpartynewie.dem,aes(x=period,y=effect,color=color,group=group)) + geom_point() + geom_errorbar(aes(ymin=lci,ymax=uci),size=0.5,width=0) + theme_classic() + theme(plot.title = element_text(hjust = 0.5)) + ggtitle("IE") + geom_hline(yintercept=0,linetype="dashed",color="red") + scale_color_identity() + xlab("Year") + ylab("Deviation from Grand Mean") + ylim(-63,28)
ggarrange(figurea50.a,figurea50.b,nrow=1)

# Figure A51
appendixcohortoutpartynewie.dem <- data.frame(cohort=c(1885,1890,1895,1900,1905,1910,1915,1920,1925,1930,1935,1940,1945,1950,1955,1960,1965,1970,1975,1980,1985,1990,1995),
                                              effect=c(-3.426,14.14,-3.222,-6.431,-0.866,0.393,0.242,0.667,0.451,1.331,2.855,0.707,-0.671,-0.847,-1.221,-1.082,0.531,0.084,-1.595,-1.787,-0.933,-1.491,2.171),
                                              se=c(9.057,6.685,3.459,2.568,1.937,1.698,1.533,1.381,1.254,1.18,1.06,0.919,0.777,0.661,0.578,0.569,0.629,0.743,0.829,0.936,1.05,1.204,1.981))
appendixcohortoutpartynewie.dem$lci <- appendixcohortoutpartynewie.dem$effect - (1.96*appendixcohortoutpartynewie.dem$se)
appendixcohortoutpartynewie.dem$uci <- appendixcohortoutpartynewie.dem$effect + (1.96*appendixcohortoutpartynewie.dem$se)
appendixcohortoutpartynewie.dem$color <- ifelse(appendixcohortoutpartynewie.dem$lci<0 & appendixcohortoutpartynewie.dem$uci>0,"Gray50","Black")
appendixcohortoutpartynewie.dem$group <- 1

appendixcohortoutpartynewapci.dem <- data.frame(cohort=c(1885,1890,1895,1900,1905,1910,1915,1920,1925,1930,1935,1940,1945,1950,1955,1960,1965,1970,1975,1980,1985,1990,1995),
                                                effect=c(as.numeric(outpartynew.dem.raw.cohort$cohort_average$cohort_average)),
                                                se=c(as.numeric(outpartynew.dem.raw.cohort$cohort_average$cohort_average_se)))
appendixcohortoutpartynewapci.dem$lci <- appendixcohortoutpartynewapci.dem$effect - (1.96*appendixcohortoutpartynewapci.dem$se)
appendixcohortoutpartynewapci.dem$uci <- appendixcohortoutpartynewapci.dem$effect + (1.96*appendixcohortoutpartynewapci.dem$se)
appendixcohortoutpartynewapci.dem$color <- ifelse(appendixcohortoutpartynewapci.dem$lci<0 & appendixcohortoutpartynewapci.dem$uci>0,"Gray50","Black")
appendixcohortoutpartynewapci.dem$group <- 1

figurea51.a <- ggplot(appendixcohortoutpartynewapci.dem,aes(x=cohort,y=effect,color=color,group=group)) + geom_point() + geom_errorbar(aes(ymin=lci,ymax=uci),size=0.5,width=0) + theme_classic() + theme(plot.title = element_text(hjust = 0.5)) + ggtitle("APCI") + geom_hline(yintercept=0,linetype="dashed",color="red") + scale_color_identity() + xlab("Birth Cohort") + ylab("Deviation from Grand Mean") + ylim(-63,28)
figurea51.b <- ggplot(appendixcohortoutpartynewie.dem,aes(x=cohort,y=effect,color=color,group=group)) + geom_point() + geom_errorbar(aes(ymin=lci,ymax=uci),size=0.5,width=0) + theme_classic() + theme(plot.title = element_text(hjust = 0.5)) + ggtitle("IE") + geom_hline(yintercept=0,linetype="dashed",color="red") + scale_color_identity() + xlab("Birth Cohort") + ylab("Deviation from Grand Mean") + ylim(-63,28)
ggarrange(figurea51.a,figurea51.b,nrow=1)

# Figure A52
appendixcohortoutpartynewapcislopes.dem <- data.frame(cohort=c(1885,1890,1895,1900,1905,1910,1915,1920,1925,1930,1935,1940,1945,1950,1955,1960,1965,1970,1975,1980,1985,1990,1995),
                                                      effect=c(as.numeric(outpartynew.dem.raw.cohort$cohort_slope$cohort_slope)),
                                                      se=c(as.numeric(outpartynew.dem.raw.cohort$cohort_slope$cohort_slope_se)))
appendixcohortoutpartynewapcislopes.dem$lci <- appendixcohortoutpartynewapcislopes.dem$effect - (1.96*appendixcohortoutpartynewapcislopes.dem$se)
appendixcohortoutpartynewapcislopes.dem$uci <- appendixcohortoutpartynewapcislopes.dem$effect + (1.96*appendixcohortoutpartynewapcislopes.dem$se)
appendixcohortoutpartynewapcislopes.dem$color <- ifelse(appendixcohortoutpartynewapcislopes.dem$lci<0 & appendixcohortoutpartynewapcislopes.dem$uci>0,"Gray50","Black")
appendixcohortoutpartynewapcislopes.dem$group <- 1

ggplot(appendixcohortoutpartynewapcislopes.dem,aes(x=cohort,y=effect,color=color,group=group)) + geom_point() + geom_errorbar(aes(ymin=lci,ymax=uci),size=0.5,width=0) + theme_classic() + theme(plot.title = element_text(hjust = 0.5)) + geom_hline(yintercept=0,linetype="dashed",color="red") + scale_color_identity() + xlab("Birth Cohort") + ylab("Deviation from Grand Mean") + ylim(-63,28)

# Figure A53
appendixageoutpartynewie.rep <- data.frame(age=c(20,25,30,35,40,45,50,55,60,65,70,75,80,85,90),
                                           effect=c(4.687,2.604,2.918,1.729,1.337,1.389,1.007,-0.536,-0.808,-2.047,-3.145,-3.666,-2.121,-2.187,-1.162),
                                           se=c(1.017,0.903,0.798,0.723,0.7,0.686,0.665,0.678,0.743,0.813,0.936,1.089,1.308,1.763,2.906))
appendixageoutpartynewie.rep$lci <- appendixageoutpartynewie.rep$effect - (1.96*appendixageoutpartynewie.rep$se)
appendixageoutpartynewie.rep$uci <- appendixageoutpartynewie.rep$effect + (1.96*appendixageoutpartynewie.rep$se)
appendixageoutpartynewie.rep$color <- ifelse(appendixageoutpartynewie.rep$lci<0 & appendixageoutpartynewie.rep$uci>0,"Gray50","Black")
appendixageoutpartynewie.rep$group <- 1

appendixageoutpartynewapci.rep <- data.frame(age=c(25,30,35,40,45,50,55,60,65,70,75,80,85),
                                             effect=c(summary(outpartynew.rep.raw$model)$coefficients[2:14,1]),
                                             se=c(summary(outpartynew.rep.raw$model)$coefficients[2:14,2]))
appendixageoutpartynewapci.rep$lci <- appendixageoutpartynewapci.rep$effect - (1.96*appendixageoutpartynewapci.rep$se)
appendixageoutpartynewapci.rep$uci <- appendixageoutpartynewapci.rep$effect + (1.96*appendixageoutpartynewapci.rep$se)
appendixageoutpartynewapci.rep$color <- ifelse(appendixageoutpartynewapci.rep$lci<0 & appendixageoutpartynewapci.rep$uci>0,"Gray50","Black")
appendixageoutpartynewapci.rep$group <- 1

figurea53.a <- ggplot(appendixageoutpartynewapci.rep,aes(x=age,y=effect,color=color,group=group)) + geom_point() + geom_errorbar(aes(ymin=lci,ymax=uci),size=0.5,width=0) + theme_classic() + theme(plot.title = element_text(hjust = 0.5)) + ggtitle("APCI") + geom_hline(yintercept=0,linetype="dashed",color="red") + scale_color_identity() + xlab("Age") + ylab("Deviation from Grand Mean") + ylim(-63,28)
figurea53.b <- ggplot(appendixageoutpartynewie.rep,aes(x=age,y=effect,color=color,group=group)) + geom_point() + geom_errorbar(aes(ymin=lci,ymax=uci),size=0.5,width=0) + theme_classic() + theme(plot.title = element_text(hjust = 0.5)) + ggtitle("IE") + geom_hline(yintercept=0,linetype="dashed",color="red") + scale_color_identity() + xlab("Age") + ylab("Deviation from Grand Mean") + ylim(-63,28)
ggarrange(figurea53.a,figurea53.b,nrow=1)

# Figure A54
appendixperiodoutpartynewie.rep <- data.frame(period=c(1975,1980,1985,1990,1995,2000,2005,2010,2015),
                                              effect=c(4.771,5.941,6.717,0.903,1.877,2.062,0.656,-9.441,-13.486),
                                              se=c(0.717,0.653,0.506,0.505,0.472,0.845,0.771,0.586,0.705))
appendixperiodoutpartynewie.rep$lci <- appendixperiodoutpartynewie.rep$effect - (1.96*appendixperiodoutpartynewie.rep$se)
appendixperiodoutpartynewie.rep$uci <- appendixperiodoutpartynewie.rep$effect + (1.96*appendixperiodoutpartynewie.rep$se)
appendixperiodoutpartynewie.rep$color <- ifelse(appendixperiodoutpartynewie.rep$lci<0 & appendixperiodoutpartynewie.rep$uci>0,"Gray50","Black")
appendixperiodoutpartynewie.rep$group <- 1

appendixperiodoutpartynewapci.rep <- data.frame(period=c(1980,1985,1990,1995,2000,2005,2010,2015),
                                                effect=c(summary(outpartynew.rep.raw$model)$coefficients[15:22,1]),
                                                se=c(summary(outpartynew.rep.raw$model)$coefficients[15:22,2]))
appendixperiodoutpartynewapci.rep$lci <- appendixperiodoutpartynewapci.rep$effect - (1.96*appendixperiodoutpartynewapci.rep$se)
appendixperiodoutpartynewapci.rep$uci <- appendixperiodoutpartynewapci.rep$effect + (1.96*appendixperiodoutpartynewapci.rep$se)
appendixperiodoutpartynewapci.rep$color <- ifelse(appendixperiodoutpartynewapci.rep$lci<0 & appendixperiodoutpartynewapci.rep$uci>0,"Gray50","Black")
appendixperiodoutpartynewapci.rep$group <- 1

figurea54.a <- ggplot(appendixperiodoutpartynewapci.rep,aes(x=period,y=effect,color=color,group=group)) + geom_point() + geom_errorbar(aes(ymin=lci,ymax=uci),size=0.5,width=0) + theme_classic() + theme(plot.title = element_text(hjust = 0.5)) + ggtitle("APCI") + geom_hline(yintercept=0,linetype="dashed",color="red") + scale_color_identity() + xlab("Year") + ylab("Deviation from Grand Mean") + ylim(-63,28)
figurea54.b <- ggplot(appendixperiodoutpartynewie.rep,aes(x=period,y=effect,color=color,group=group)) + geom_point() + geom_errorbar(aes(ymin=lci,ymax=uci),size=0.5,width=0) + theme_classic() + theme(plot.title = element_text(hjust = 0.5)) + ggtitle("IE") + geom_hline(yintercept=0,linetype="dashed",color="red") + scale_color_identity() + xlab("Year") + ylab("Deviation from Grand Mean") + ylim(-63,28)
ggarrange(figurea54.a,figurea54.b,nrow=1)

# Figure A55
appendixcohortoutpartynewie.rep <- data.frame(cohort=c(1885,1890,1895,1900,1905,1910,1915,1920,1925,1930,1935,1940,1945,1950,1955,1960,1965,1970,1975,1980,1985,1990,1995),
                                              effect=c(-42.671,16.678,1.332,7.563,4.969,6.336,4.991,6.274,3.862,-0.256,0.576,-0.575,-1.775,-1.819,-3.788,-1.725,-1.991,-0.124,-1.502,-0.684,0.885,0.106,3.338),
                                              se=c(10.281,5.215,3.102,2.61,2.08,1.785,1.636,1.486,1.302,1.208,1.079,0.938,0.811,0.706,0.607,0.582,0.65,0.784,0.95,1.145,1.176,1.363,2.132))
appendixcohortoutpartynewie.rep$lci <- appendixcohortoutpartynewie.rep$effect - (1.96*appendixcohortoutpartynewie.rep$se)
appendixcohortoutpartynewie.rep$uci <- appendixcohortoutpartynewie.rep$effect + (1.96*appendixcohortoutpartynewie.rep$se)
appendixcohortoutpartynewie.rep$color <- ifelse(appendixcohortoutpartynewie.rep$lci<0 & appendixcohortoutpartynewie.rep$uci>0,"Gray50","Black")
appendixcohortoutpartynewie.rep$group <- 1

appendixcohortoutpartynewapci.rep <- data.frame(cohort=c(1890,1895,1900,1905,1910,1915,1920,1925,1930,1935,1940,1945,1950,1955,1960,1965,1970,1975,1980,1985,1990,1995),
                                                effect=c(as.numeric(outpartynew.rep.raw.cohort$cohort_average$cohort_average)),
                                                se=c(as.numeric(outpartynew.rep.raw.cohort$cohort_average$cohort_average_se)))
appendixcohortoutpartynewapci.rep$lci <- appendixcohortoutpartynewapci.rep$effect - (1.96*appendixcohortoutpartynewapci.rep$se)
appendixcohortoutpartynewapci.rep$uci <- appendixcohortoutpartynewapci.rep$effect + (1.96*appendixcohortoutpartynewapci.rep$se)
appendixcohortoutpartynewapci.rep$color <- ifelse(appendixcohortoutpartynewapci.rep$lci<0 & appendixcohortoutpartynewapci.rep$uci>0,"Gray50","Black")
appendixcohortoutpartynewapci.rep$group <- 1

figurea55.a <- ggplot(appendixcohortoutpartynewapci.rep,aes(x=cohort,y=effect,color=color,group=group)) + geom_point() + geom_errorbar(aes(ymin=lci,ymax=uci),size=0.5,width=0) + theme_classic() + theme(plot.title = element_text(hjust = 0.5)) + ggtitle("APCI") + geom_hline(yintercept=0,linetype="dashed",color="red") + scale_color_identity() + xlab("Birth Cohort") + ylab("Deviation from Grand Mean") + ylim(-63,28)
figurea55.b <- ggplot(appendixcohortoutpartynewie.rep,aes(x=cohort,y=effect,color=color,group=group)) + geom_point() + geom_errorbar(aes(ymin=lci,ymax=uci),size=0.5,width=0) + theme_classic() + theme(plot.title = element_text(hjust = 0.5)) + ggtitle("IE") + geom_hline(yintercept=0,linetype="dashed",color="red") + scale_color_identity() + xlab("Birth Cohort") + ylab("Deviation from Grand Mean") + ylim(-63,28)
ggarrange(figurea55.a,figurea55.b,nrow=1)

# Figure A56
appendixcohortoutpartynewapcislopes.rep <- data.frame(cohort=c(1890,1895,1900,1905,1910,1915,1920,1925,1930,1935,1940,1945,1950,1955,1960,1965,1970,1975,1980,1985,1990,1995),
                                                      effect=c(as.numeric(outpartynew.rep.raw.cohort$cohort_slope$cohort_slope)),
                                                      se=c(as.numeric(outpartynew.rep.raw.cohort$cohort_slope$cohort_slope_se)))
appendixcohortoutpartynewapcislopes.rep$lci <- appendixcohortoutpartynewapcislopes.rep$effect - (1.96*appendixcohortoutpartynewapcislopes.rep$se)
appendixcohortoutpartynewapcislopes.rep$uci <- appendixcohortoutpartynewapcislopes.rep$effect + (1.96*appendixcohortoutpartynewapcislopes.rep$se)
appendixcohortoutpartynewapcislopes.rep$color <- ifelse(appendixcohortoutpartynewapcislopes.rep$lci<0 & appendixcohortoutpartynewapcislopes.rep$uci>0,"Gray50","Black")
appendixcohortoutpartynewapcislopes.rep$group <- 1

ggplot(appendixcohortoutpartynewapcislopes.rep,aes(x=cohort,y=effect,color=color,group=group)) + geom_point() + geom_errorbar(aes(ymin=lci,ymax=uci),size=0.5,width=0) + theme_classic() + theme(plot.title = element_text(hjust = 0.5)) + geom_hline(yintercept=0,linetype="dashed",color="red") + scale_color_identity() + xlab("Birth Cohort") + ylab("Deviation from Grand Mean") + ylim(-63,28)

# Table A33
## Model 1
outpartycombined.all.raw <- APCI::temp_model(data=anes.cumulative.full,outcome='outparty.combined',age='agegroup',period='yeargroup.old',cohort='cohort.combined',covariate=c(),family='gaussian')
outpartycombined.all.raw.cohort <- APCI::cohortdeviation(outpartycombined.all.raw$A,outpartycombined.all.raw$P,outpartycombined.all.raw$C,model=outpartycombined.all.raw$model,covariate=c())
## Model 2
outpartycombinedpid.all.raw <- APCI::temp_model(data=anes.cumulative.full,outcome='outparty.combined',age='agegroup',period='yeargroup.old',cohort='cohort.combined',covariate=c('strong'),family='gaussian')
outpartycombinedpid.all.raw.cohort <- APCI::cohortdeviation(outpartycombinedpid.all.raw$A,outpartycombinedpid.all.raw$P,outpartycombinedpid.all.raw$C,model=outpartycombinedpid.all.raw$model,covariate=c('strong'))
## Model 3
outpartycombinedsort.all.raw <- APCI::temp_model(data=repnot90,outcome='outparty.combined',age='agegroup',period='yeargroup.combined',cohort='cohort.combined',covariate=c('ideosort','socialsort'),family='gaussian')
outpartycombinedsort.all.raw.cohort <- APCI::cohortdeviation(outpartycombinedsort.all.raw$A,outpartycombinedsort.all.raw$P,outpartycombinedsort.all.raw$C,model=outpartycombinedsort.all.raw$model,covariate=c('strong'))
## Model 4
outpartycombined.all.controls <- APCI::temp_model(data=repnot90,outcome='outparty.combined',age='agegroup',period='yeargroup.combined',cohort='cohort.combined',covariate=c('gender','black','education','ideo','attendance','south'),family='gaussian')
outpartycombined.all.controls.cohort <- APCI::cohortdeviation(outpartycombined.all.controls$A,outpartycombined.all.controls$P,outpartycombined.all.controls$C,model=outpartycombined.all.controls$model,covariate=c('gender','black','education','ideo','attendance','south'))
## Model 5
outpartycombinedpid.all.controls <- APCI::temp_model(data=repnot90,outcome='outparty.combined',age='agegroup',period='yeargroup.combined',cohort='cohort.combined',covariate=c('strong','gender','black','education','ideo','attendance','south'),family='gaussian')
outpartycombinedpid.all.controls.cohort <- APCI::cohortdeviation(outpartycombinedpid.all.controls$A,outpartycombinedpid.all.controls$P,outpartycombinedpid.all.controls$C,model=outpartycombinedpid.all.controls$model,covariate=c('strong','gender','black','education','ideo','attendance','south'))
## Model 6
outpartycombinedsort.all.controls <- APCI::temp_model(data=repnot90,outcome='outparty.combined',age='agegroup',period='yeargroup.combined',cohort='cohort.combined',covariate=c('ideosort','socialsort','gender','black','education','ideo','attendance','south'),family='gaussian')
outpartycombinedsort.all.controls.cohort <- APCI::cohortdeviation(outpartycombinedsort.all.controls$A,outpartycombinedsort.all.controls$P,outpartycombinedsort.all.controls$C,model=outpartycombinedsort.all.controls$model,covariate=c('ideosort','socialsort','gender','black','education','ideo','attendance','south'))

# Table A35
## Model 1
outpartycombined.dem.raw <- APCI::temp_model(data=anes.cumulative.dem,outcome='outparty.combined',age='agegroup',period='yeargroup.old',cohort='cohort.combined',covariate=c(),family='gaussian')
outpartycombined.dem.raw.cohort <- APCI::cohortdeviation(outpartycombined.dem.raw$A,outpartycombined.dem.raw$P,outpartycombined.dem.raw$C,model=outpartycombined.dem.raw$model,covariate=c())
## Model 2
outpartycombinedpid.dem.raw <- APCI::temp_model(data=anes.cumulative.dem,outcome='outparty.combined',age='agegroup',period='yeargroup.old',cohort='cohort.combined',covariate=c('strong'),family='gaussian')
outpartycombinedpid.dem.raw.cohort <- APCI::cohortdeviation(outpartycombinedpid.dem.raw$A,outpartycombinedpid.dem.raw$P,outpartycombinedpid.dem.raw$C,model=outpartycombinedpid.dem.raw$model,covariate=c('strong'))
## Model 3
outpartycombinedsort.dem.raw <- APCI::temp_model(data=demnot90,outcome='outparty.combined',age='agegroup',period='yeargroup.combined',cohort='cohort.combined',covariate=c('ideosort','socialsort'),family='gaussian')
outpartycombinedsort.dem.raw.cohort <- APCI::cohortdeviation(outpartycombinedsort.dem.raw$A,outpartycombinedsort.dem.raw$P,outpartycombinedsort.dem.raw$C,model=outpartycombinedsort.dem.raw$model,covariate=c('strong'))
## Model 4
outpartycombined.dem.controls <- APCI::temp_model(data=demnot90,outcome='outparty.combined',age='agegroup',period='yeargroup.combined',cohort='cohort.combined',covariate=c('gender','black','education','ideo','attendance','south'),family='gaussian')
outpartycombined.dem.controls.cohort <- APCI::cohortdeviation(outpartycombined.dem.controls$A,outpartycombined.dem.controls$P,outpartycombined.dem.controls$C,model=outpartycombined.dem.controls$model,covariate=c('gender','black','education','ideo','attendance','south'))
## Model 5
outpartycombinedpid.dem.controls <- APCI::temp_model(data=demnot90,outcome='outparty.combined',age='agegroup',period='yeargroup.combined',cohort='cohort.combined',covariate=c('strong','gender','black','education','ideo','attendance','south'),family='gaussian')
outpartycombinedpid.dem.controls.cohort <- APCI::cohortdeviation(outpartycombinedpid.dem.controls$A,outpartycombinedpid.dem.controls$P,outpartycombinedpid.dem.controls$C,model=outpartycombinedpid.dem.controls$model,covariate=c('strong','gender','black','education','ideo','attendance','south'))
## Model 6
outpartycombinedsort.dem.controls <- APCI::temp_model(data=demnot90,outcome='outparty.combined',age='agegroup',period='yeargroup.combined',cohort='cohort.combined',covariate=c('ideosort','socialsort','gender','black','education','ideo','attendance','south'),family='gaussian')
outpartycombinedsort.dem.controls.cohort <- APCI::cohortdeviation(outpartycombinedsort.dem.controls$A,outpartycombinedsort.dem.controls$P,outpartycombinedsort.dem.controls$C,model=outpartycombinedsort.dem.controls$model,covariate=c('ideosort','socialsort','gender','black','education','ideo','attendance','south'))

# Table A37
## Model 1
outpartycombined.rep.raw <- APCI::temp_model(data=repnot90,outcome='outparty.combined',age='agegroup',period='yeargroup.old',cohort='cohort.combined',covariate=c(),family='gaussian')
outpartycombined.rep.raw.cohort <- APCI::cohortdeviation(outpartycombined.rep.raw$A,outpartycombined.rep.raw$P,outpartycombined.rep.raw$C,model=outpartycombined.rep.raw$model,covariate=c())
## Model 2
outpartycombinedpid.rep.raw <- APCI::temp_model(data=repnot90,outcome='outparty.combined',age='agegroup',period='yeargroup.old',cohort='cohort.combined',covariate=c('strong'),family='gaussian')
outpartycombinedpid.rep.raw.cohort <- APCI::cohortdeviation(outpartycombinedpid.rep.raw$A,outpartycombinedpid.rep.raw$P,outpartycombinedpid.rep.raw$C,model=outpartycombinedpid.rep.raw$model,covariate=c('strong'))
## Model 3
outpartycombinedsort.rep.raw <- APCI::temp_model(data=repnot90,outcome='outparty.combined',age='agegroup',period='yeargroup.combined',cohort='cohort.combined',covariate=c('ideosort','socialsort'),family='gaussian')
outpartycombinedsort.rep.raw.cohort <- APCI::cohortdeviation(outpartycombinedsort.rep.raw$A,outpartycombinedsort.rep.raw$P,outpartycombinedsort.rep.raw$C,model=outpartycombinedsort.rep.raw$model,covariate=c('strong'))
## Model 4
outpartycombined.rep.controls <- APCI::temp_model(data=repnot90,outcome='outparty.combined',age='agegroup',period='yeargroup.combined',cohort='cohort.combined',covariate=c('gender','black','education','ideo','attendance','south'),family='gaussian')
outpartycombined.rep.controls.cohort <- APCI::cohortdeviation(outpartycombined.rep.controls$A,outpartycombined.rep.controls$P,outpartycombined.rep.controls$C,model=outpartycombined.rep.controls$model,covariate=c('gender','black','education','ideo','attendance','south'))
## Model 5
outpartycombinedpid.rep.controls <- APCI::temp_model(data=repnot90,outcome='outparty.combined',age='agegroup',period='yeargroup.combined',cohort='cohort.combined',covariate=c('strong','gender','black','education','ideo','attendance','south'),family='gaussian')
outpartycombinedpid.rep.controls.cohort <- APCI::cohortdeviation(outpartycombinedpid.rep.controls$A,outpartycombinedpid.rep.controls$P,outpartycombinedpid.rep.controls$C,model=outpartycombinedpid.rep.controls$model,covariate=c('strong','gender','black','education','ideo','attendance','south'))
## Model 6
outpartycombinedsort.rep.controls <- APCI::temp_model(data=repnot90,outcome='outparty.combined',age='agegroup',period='yeargroup.combined',cohort='cohort.combined',covariate=c('ideosort','socialsort','gender','black','education','ideo','attendance','south'),family='gaussian')
outpartycombinedsort.rep.controls.cohort <- APCI::cohortdeviation(outpartycombinedsort.rep.controls$A,outpartycombinedsort.rep.controls$P,outpartycombinedsort.rep.controls$C,model=outpartycombinedsort.rep.controls$model,covariate=c('ideosort','socialsort','gender','black','education','ideo','attendance','south'))

# Figure A57
appendixageoutpartycombinedie.all <- data.frame(age=c(20,25,30,35,40,45,50,55,60,65,70,75,80,85,90),
                                                effect=c(1.459,0.786,1.411,0.662,0.276,-0.035,0.854,0.113,0.651,-0.271,-0.966,0.361,-0.158,-1.505,-3.637),
                                                se=c(0.742,0.65,0.571,0.501,0.454,0.422,0.401,0.409,0.445,0.5,0.592,0.705,0.865,1.202,2.007))
appendixageoutpartycombinedie.all$lci <- appendixageoutpartycombinedie.all$effect - (1.96*appendixageoutpartycombinedie.all$se)
appendixageoutpartycombinedie.all$uci <- appendixageoutpartycombinedie.all$effect + (1.96*appendixageoutpartycombinedie.all$se)
appendixageoutpartycombinedie.all$color <- ifelse(appendixageoutpartycombinedie.all$lci<0 & appendixageoutpartycombinedie.all$uci>0,"Gray50","Black")
appendixageoutpartycombinedie.all$group <- 1

appendixageoutpartycombinedapci.all <- data.frame(age=c(25,30,35,40,45,50,55,60,65,70,75,80,85,90),
                                                  effect=c(summary(outpartycombined.all.raw$model)$coefficients[2:15,1]),
                                                  se=c(summary(outpartycombined.all.raw$model)$coefficients[2:15,2]))
appendixageoutpartycombinedapci.all$lci <- appendixageoutpartycombinedapci.all$effect - (1.96*appendixageoutpartycombinedapci.all$se)
appendixageoutpartycombinedapci.all$uci <- appendixageoutpartycombinedapci.all$effect + (1.96*appendixageoutpartycombinedapci.all$se)
appendixageoutpartycombinedapci.all$color <- ifelse(appendixageoutpartycombinedapci.all$lci<0 & appendixageoutpartycombinedapci.all$uci>0,"Gray50","Black")
appendixageoutpartycombinedapci.all$group <- 1

figurea57.a <- ggplot(appendixageoutpartycombinedapci.all,aes(x=age,y=effect,color=color,group=group)) + geom_point() + geom_errorbar(aes(ymin=lci,ymax=uci),size=0.5,width=0) + theme_classic() + theme(plot.title = element_text(hjust = 0.5)) + ggtitle("APCI") + geom_hline(yintercept=0,linetype="dashed",color="red") + scale_color_identity() + xlab("Age") + ylab("Deviation from Grand Mean") + ylim(-31,41)
figurea57.b <- ggplot(appendixageoutpartycombinedie.all,aes(x=age,y=effect,color=color,group=group)) + geom_point() + geom_errorbar(aes(ymin=lci,ymax=uci),size=0.5,width=0) + theme_classic() + theme(plot.title = element_text(hjust = 0.5)) + ggtitle("IE") + geom_hline(yintercept=0,linetype="dashed",color="red") + scale_color_identity() + xlab("Age") + ylab("Deviation from Grand Mean") + ylim(-31,41)
ggarrange(figurea57.a,figurea57.b,nrow=1)

# Figure A58
appendixperiodoutpartycombinedie.all <- data.frame(period=c(1960,1965,1970,1975,1980,1985,1990,1995,2000,2005,2010,2015),
                                                   effect=c(7.039,9.996,12.742,6.538,2.875,3.866,-0.481,-1.257,-3.291,-6.789,-14.756,-16.481),
                                                   se=c(0.716,0.554,0.509,0.403,0.404,0.306,0.348,0.355,0.65,0.564,0.532,0.648))
appendixperiodoutpartycombinedie.all$lci <- appendixperiodoutpartycombinedie.all$effect - (1.96*appendixperiodoutpartycombinedie.all$se)
appendixperiodoutpartycombinedie.all$uci <- appendixperiodoutpartycombinedie.all$effect + (1.96*appendixperiodoutpartycombinedie.all$se)
appendixperiodoutpartycombinedie.all$color <- ifelse(appendixperiodoutpartycombinedie.all$lci<0 & appendixperiodoutpartycombinedie.all$uci>0,"Gray50","Black")
appendixperiodoutpartycombinedie.all$group <- 1

appendixperiodoutpartycombinedapci.all <- data.frame(period=c(1965,1970,1975,1980,1985,1990,1995,2000,2005,2010,2015),
                                                     effect=c(summary(outpartycombined.all.raw$model)$coefficients[16:26,1]),
                                                     se=c(summary(outpartycombined.all.raw$model)$coefficients[16:26,2]))
appendixperiodoutpartycombinedapci.all$lci <- appendixperiodoutpartycombinedapci.all$effect - (1.96*appendixperiodoutpartycombinedapci.all$se)
appendixperiodoutpartycombinedapci.all$uci <- appendixperiodoutpartycombinedapci.all$effect + (1.96*appendixperiodoutpartycombinedapci.all$se)
appendixperiodoutpartycombinedapci.all$color <- ifelse(appendixperiodoutpartycombinedapci.all$lci<0 & appendixperiodoutpartycombinedapci.all$uci>0,"Gray50","Black")
appendixperiodoutpartycombinedapci.all$group <- 1

figurea58.a <- ggplot(appendixperiodoutpartycombinedapci.all,aes(x=period,y=effect,color=color,group=group)) + geom_point() + geom_errorbar(aes(ymin=lci,ymax=uci),size=0.5,width=0) + theme_classic() + theme(plot.title = element_text(hjust = 0.5)) + ggtitle("APCI") + geom_hline(yintercept=0,linetype="dashed",color="red") + scale_color_identity() + xlab("Year") + ylab("Deviation from Grand Mean") + ylim(-31,41)
figurea58.b <- ggplot(appendixperiodoutpartycombinedie.all,aes(x=period,y=effect,color=color,group=group)) + geom_point() + geom_errorbar(aes(ymin=lci,ymax=uci),size=0.5,width=0) + theme_classic() + theme(plot.title = element_text(hjust = 0.5)) + ggtitle("IE") + geom_hline(yintercept=0,linetype="dashed",color="red") + scale_color_identity() + xlab("Year") + ylab("Deviation from Grand Mean") + ylim(-31,41)
ggarrange(figurea58.a,figurea58.b,nrow=1)

# Figure A59
appendixcohortoutpartycombinedie.all <- data.frame(cohort=c(1870,1875,1880,1885,1890,1895,1900,1905,1910,1915,1920,1925,1930,1935,1940,1945,1950,1955,1960,1965,1970,1975,1980,1985,1990,1995),
                                                   effect=c(2.681,1.488,1.207,2.802,3.53,0.396,-0.034,0.537,0.974,0.674,0.166,-0.407,-0.95,-0.267,-1.779,-2.508,-2.77,-3.262,-1.981,-0.945,-0.182,-1.655,-1.416,0.098,-0.344,3.946),
                                                   se=c(9.382,6.395,3.613,2.208,1.81,1.544,1.398,1.245,1.122,1.018,0.918,0.82,0.744,0.657,0.562,0.477,0.414,0.381,0.393,0.459,0.562,0.662,0.77,0.849,0.975,1.504))
appendixcohortoutpartycombinedie.all$lci <- appendixcohortoutpartycombinedie.all$effect - (1.96*appendixcohortoutpartycombinedie.all$se)
appendixcohortoutpartycombinedie.all$uci <- appendixcohortoutpartycombinedie.all$effect + (1.96*appendixcohortoutpartycombinedie.all$se)
appendixcohortoutpartycombinedie.all$color <- ifelse(appendixcohortoutpartycombinedie.all$lci<0 & appendixcohortoutpartycombinedie.all$uci>0,"Gray50","Black")
appendixcohortoutpartycombinedie.all$group <- 1

appendixcohortoutpartycombinedapci.all <- data.frame(cohort=c(1870,1875,1880,1885,1890,1895,1900,1905,1910,1915,1920,1925,1930,1935,1940,1945,1950,1955,1960,1965,1970,1975,1980,1985,1990,1995),
                                                     effect=c(as.numeric(outpartycombined.all.raw.cohort$cohort_average$cohort_average)),
                                                     se=c(as.numeric(outpartycombined.all.raw.cohort$cohort_average$cohort_average_se)))
appendixcohortoutpartycombinedapci.all$lci <- appendixcohortoutpartycombinedapci.all$effect - (1.96*appendixcohortoutpartycombinedapci.all$se)
appendixcohortoutpartycombinedapci.all$uci <- appendixcohortoutpartycombinedapci.all$effect + (1.96*appendixcohortoutpartycombinedapci.all$se)
appendixcohortoutpartycombinedapci.all$color <- ifelse(appendixcohortoutpartycombinedapci.all$lci<0 & appendixcohortoutpartycombinedapci.all$uci>0,"Gray50","Black")
appendixcohortoutpartycombinedapci.all$group <- 1

figurea59.a <- ggplot(appendixcohortoutpartycombinedapci.all,aes(x=cohort,y=effect,color=color,group=group)) + geom_point() + geom_errorbar(aes(ymin=lci,ymax=uci),size=0.5,width=0) + theme_classic() + theme(plot.title = element_text(hjust = 0.5)) + ggtitle("APCI") + geom_hline(yintercept=0,linetype="dashed",color="red") + scale_color_identity() + xlab("Birth Cohort") + ylab("Deviation from Grand Mean") + ylim(-31,41)
figurea59.b <- ggplot(appendixcohortoutpartycombinedie.all,aes(x=cohort,y=effect,color=color,group=group)) + geom_point() + geom_errorbar(aes(ymin=lci,ymax=uci),size=0.5,width=0) + theme_classic() + theme(plot.title = element_text(hjust = 0.5)) + ggtitle("IE") + geom_hline(yintercept=0,linetype="dashed",color="red") + scale_color_identity() + xlab("Birth Cohort") + ylab("Deviation from Grand Mean") + ylim(-31,41)
ggarrange(figurea59.a,figurea59.b,nrow=1)

# Figure A60
appendixcohortoutpartycombinedapcislopes.all <- data.frame(cohort=c(1870,1875,1880,1885,1890,1895,1900,1905,1910,1915,1920,1925,1930,1935,1940,1945,1950,1955,1960,1965,1970,1975,1980,1985,1990,1995),
                                                           effect=c(as.numeric(outpartycombined.all.raw.cohort$cohort_slope$cohort_slope)),
                                                           se=c(as.numeric(outpartycombined.all.raw.cohort$cohort_slope$cohort_slope_se)))
appendixcohortoutpartycombinedapcislopes.all$lci <- appendixcohortoutpartycombinedapcislopes.all$effect - (1.96*appendixcohortoutpartycombinedapcislopes.all$se)
appendixcohortoutpartycombinedapcislopes.all$uci <- appendixcohortoutpartycombinedapcislopes.all$effect + (1.96*appendixcohortoutpartycombinedapcislopes.all$se)
appendixcohortoutpartycombinedapcislopes.all$color <- ifelse(appendixcohortoutpartycombinedapcislopes.all$lci<0 & appendixcohortoutpartycombinedapcislopes.all$uci>0,"Gray50","Black")
appendixcohortoutpartycombinedapcislopes.all$group <- 1

ggplot(appendixcohortoutpartycombinedapcislopes.all,aes(x=cohort,y=effect,color=color,group=group)) + geom_point() + geom_errorbar(aes(ymin=lci,ymax=uci),size=0.5,width=0) + theme_classic() + theme(plot.title = element_text(hjust = 0.5)) + geom_hline(yintercept=0,linetype="dashed",color="red") + scale_color_identity() + xlab("Birth Cohort") + ylab("Deviation from Grand Mean") + ylim(-31,41)

# Figure A61
appendixageoutpartycombinedie.dem <- data.frame(age=c(20,25,30,35,40,45,50,55,60,65,70,75,80,85,90),
                                                effect=c(1.456,1.121,1.871,0.974,0.654,-0.298,1.110,0.489,0.817,0.093,-0.765,1.489,0.155,-2.913,-6.253),
                                                se=c(1.332,1.157,0.998,0.846,0.725,0.627,0.563,0.561,0.627,0.73,0.9,1.101,1.357,1.825,3.135))
appendixageoutpartycombinedie.dem$lci <- appendixageoutpartycombinedie.dem$effect - (1.96*appendixageoutpartycombinedie.dem$se)
appendixageoutpartycombinedie.dem$uci <- appendixageoutpartycombinedie.dem$effect + (1.96*appendixageoutpartycombinedie.dem$se)
appendixageoutpartycombinedie.dem$color <- ifelse(appendixageoutpartycombinedie.dem$lci<0 & appendixageoutpartycombinedie.dem$uci>0,"Gray50","Black")
appendixageoutpartycombinedie.dem$group <- 1

appendixageoutpartycombinedapci.dem <- data.frame(age=c(25,30,35,40,45,50,55,60,65,70,75,80,85,90),
                                                  effect=c(summary(outpartycombined.dem.raw$model)$coefficients[2:15,1]),
                                                  se=c(summary(outpartycombined.dem.raw$model)$coefficients[2:15,2]))
appendixageoutpartycombinedapci.dem$lci <- appendixageoutpartycombinedapci.dem$effect - (1.96*appendixageoutpartycombinedapci.dem$se)
appendixageoutpartycombinedapci.dem$uci <- appendixageoutpartycombinedapci.dem$effect + (1.96*appendixageoutpartycombinedapci.dem$se)
appendixageoutpartycombinedapci.dem$color <- ifelse(appendixageoutpartycombinedapci.dem$lci<0 & appendixageoutpartycombinedapci.dem$uci>0,"Gray50","Black")
appendixageoutpartycombinedapci.dem$group <- 1

figurea61.a <- ggplot(appendixageoutpartycombinedapci.dem,aes(x=age,y=effect,color=color,group=group)) + geom_point() + geom_errorbar(aes(ymin=lci,ymax=uci),size=0.5,width=0) + theme_classic() + theme(plot.title = element_text(hjust = 0.5)) + ggtitle("APCI") + geom_hline(yintercept=0,linetype="dashed",color="red") + scale_color_identity() + xlab("Age") + ylab("Deviation from Grand Mean") + ylim(-31,41)
figurea61.b <- ggplot(appendixageoutpartycombinedie.dem,aes(x=age,y=effect,color=color,group=group)) + geom_point() + geom_errorbar(aes(ymin=lci,ymax=uci),size=0.5,width=0) + theme_classic() + theme(plot.title = element_text(hjust = 0.5)) + ggtitle("IE") + geom_hline(yintercept=0,linetype="dashed",color="red") + scale_color_identity() + xlab("Age") + ylab("Deviation from Grand Mean") + ylim(-31,41)
ggarrange(figurea61.a,figurea61.b,nrow=1)

# Figure A62
appendixperiodoutpartycombinedie.dem <- data.frame(period=c(1960,1965,1970,1975,1980,1985,1990,1995,2000,2005,2010,2015),
                                                   effect=c(6.525,11.096,12.398,7.128,2.045,3.537,0.721,-1.158,-4.513,-8.094,-15.017,-14.670),
                                                   se=c(1.171,0.934,0.816,0.625,0.579,0.422,0.488,0.52,0.959,0.857,0.917,1.133))
appendixperiodoutpartycombinedie.dem$lci <- appendixperiodoutpartycombinedie.dem$effect - (1.96*appendixperiodoutpartycombinedie.dem$se)
appendixperiodoutpartycombinedie.dem$uci <- appendixperiodoutpartycombinedie.dem$effect + (1.96*appendixperiodoutpartycombinedie.dem$se)
appendixperiodoutpartycombinedie.dem$color <- ifelse(appendixperiodoutpartycombinedie.dem$lci<0 & appendixperiodoutpartycombinedie.dem$uci>0,"Gray50","Black")
appendixperiodoutpartycombinedie.dem$group <- 1

appendixperiodoutpartycombinedapci.dem <- data.frame(period=c(1965,1970,1975,1980,1985,1990,1995,2000,2005,2010,2015),
                                                     effect=c(summary(outpartycombined.dem.raw$model)$coefficients[16:26,1]),
                                                     se=c(summary(outpartycombined.dem.raw$model)$coefficients[16:26,2]))
appendixperiodoutpartycombinedapci.dem$lci <- appendixperiodoutpartycombinedapci.dem$effect - (1.96*appendixperiodoutpartycombinedapci.dem$se)
appendixperiodoutpartycombinedapci.dem$uci <- appendixperiodoutpartycombinedapci.dem$effect + (1.96*appendixperiodoutpartycombinedapci.dem$se)
appendixperiodoutpartycombinedapci.dem$color <- ifelse(appendixperiodoutpartycombinedapci.dem$lci<0 & appendixperiodoutpartycombinedapci.dem$uci>0,"Gray50","Black")
appendixperiodoutpartycombinedapci.dem$group <- 1

figurea62.a <- ggplot(appendixperiodoutpartycombinedapci.dem,aes(x=period,y=effect,color=color,group=group)) + geom_point() + geom_errorbar(aes(ymin=lci,ymax=uci),size=0.5,width=0) + theme_classic() + theme(plot.title = element_text(hjust = 0.5)) + ggtitle("APCI") + geom_hline(yintercept=0,linetype="dashed",color="red") + scale_color_identity() + xlab("Year") + ylab("Deviation from Grand Mean") + ylim(-31,41)
figurea62.b <- ggplot(appendixperiodoutpartycombinedie.dem,aes(x=period,y=effect,color=color,group=group)) + geom_point() + geom_errorbar(aes(ymin=lci,ymax=uci),size=0.5,width=0) + theme_classic() + theme(plot.title = element_text(hjust = 0.5)) + ggtitle("IE") + geom_hline(yintercept=0,linetype="dashed",color="red") + scale_color_identity() + xlab("Year") + ylab("Deviation from Grand Mean") + ylim(-31,41)
ggarrange(figurea62.a,figurea62.b,nrow=1)

# Figure A63
appendixcohortoutpartycombinedie.dem <- data.frame(cohort=c(1870,1875,1880,1885,1890,1895,1900,1905,1910,1915,1920,1925,1930,1935,1940,1945,1950,1955,1960,1965,1970,1975,1980,1985,1990,1995),
                                                   effect=c(7.731,-4.667,-1.586,5.148,1.794,1.889,-1.216,1.45,0.65,0.361,0.35,0.096,0.456,1.301,-0.633,-1.769,-2.19,-1.936,-1.704,-0.062,0.602,-2.068,-2.155,-1.261,-1.647,2.268),
                                                   se=c(1.231,10.075,5.533,3.596,3.069,2.663,2.409,2.162,1.955,1.768,1.579,1.403,1.244,1.082,0.904,0.737,0.601,0.53,0.545,0.651,0.817,0.965,1.131,1.308,1.511,2.226))
appendixcohortoutpartycombinedie.dem$lci <- appendixcohortoutpartycombinedie.dem$effect - (1.96*appendixcohortoutpartycombinedie.dem$se)
appendixcohortoutpartycombinedie.dem$uci <- appendixcohortoutpartycombinedie.dem$effect + (1.96*appendixcohortoutpartycombinedie.dem$se)
appendixcohortoutpartycombinedie.dem$color <- ifelse(appendixcohortoutpartycombinedie.dem$lci<0 & appendixcohortoutpartycombinedie.dem$uci>0,"Gray50","Black")
appendixcohortoutpartycombinedie.dem$group <- 1

appendixcohortoutpartycombinedapci.dem <- data.frame(cohort=c(1870,1875,1880,1885,1890,1895,1900,1905,1910,1915,1920,1925,1930,1935,1940,1945,1950,1955,1960,1965,1970,1975,1980,1985,1990,1995),
                                                     effect=c(as.numeric(outpartycombined.dem.raw.cohort$cohort_average$cohort_average)),
                                                     se=c(as.numeric(outpartycombined.dem.raw.cohort$cohort_average$cohort_average_se)))
appendixcohortoutpartycombinedapci.dem$lci <- appendixcohortoutpartycombinedapci.dem$effect - (1.96*appendixcohortoutpartycombinedapci.dem$se)
appendixcohortoutpartycombinedapci.dem$uci <- appendixcohortoutpartycombinedapci.dem$effect + (1.96*appendixcohortoutpartycombinedapci.dem$se)
appendixcohortoutpartycombinedapci.dem$color <- ifelse(appendixcohortoutpartycombinedapci.dem$lci<0 & appendixcohortoutpartycombinedapci.dem$uci>0,"Gray50","Black")
appendixcohortoutpartycombinedapci.dem$group <- 1

figurea63.a <- ggplot(appendixcohortoutpartycombinedapci.dem,aes(x=cohort,y=effect,color=color,group=group)) + geom_point() + geom_errorbar(aes(ymin=lci,ymax=uci),size=0.5,width=0) + theme_classic() + theme(plot.title = element_text(hjust = 0.5)) + ggtitle("APCI") + geom_hline(yintercept=0,linetype="dashed",color="red") + scale_color_identity() + xlab("Birth Cohort") + ylab("Deviation from Grand Mean") + ylim(-31,41)
figurea63.b <- ggplot(appendixcohortoutpartycombinedie.dem,aes(x=cohort,y=effect,color=color,group=group)) + geom_point() + geom_errorbar(aes(ymin=lci,ymax=uci),size=0.5,width=0) + theme_classic() + theme(plot.title = element_text(hjust = 0.5)) + ggtitle("IE") + geom_hline(yintercept=0,linetype="dashed",color="red") + scale_color_identity() + xlab("Birth Cohort") + ylab("Deviation from Grand Mean") + ylim(-31,41)
ggarrange(figurea63.a,figurea63.b,nrow=1)

# Figure A64
appendixcohortoutpartycombinedapcislopes.dem <- data.frame(cohort=c(1870,1875,1880,1885,1890,1895,1900,1905,1910,1915,1920,1925,1930,1935,1940,1945,1950,1955,1960,1965,1970,1975,1980,1985,1990,1995),
                                                           effect=c(as.numeric(outpartycombined.dem.raw.cohort$cohort_slope$cohort_slope)),
                                                           se=c(as.numeric(outpartycombined.dem.raw.cohort$cohort_slope$cohort_slope_se)))
appendixcohortoutpartycombinedapcislopes.dem$lci <- appendixcohortoutpartycombinedapcislopes.dem$effect - (1.96*appendixcohortoutpartycombinedapcislopes.dem$se)
appendixcohortoutpartycombinedapcislopes.dem$uci <- appendixcohortoutpartycombinedapcislopes.dem$effect + (1.96*appendixcohortoutpartycombinedapcislopes.dem$se)
appendixcohortoutpartycombinedapcislopes.dem$color <- ifelse(appendixcohortoutpartycombinedapcislopes.dem$lci<0 & appendixcohortoutpartycombinedapcislopes.dem$uci>0,"Gray50","Black")
appendixcohortoutpartycombinedapcislopes.dem$group <- 1

ggplot(appendixcohortoutpartycombinedapcislopes.dem,aes(x=cohort,y=effect,color=color,group=group)) + geom_point() + geom_errorbar(aes(ymin=lci,ymax=uci),size=0.5,width=0) + theme_classic() + theme(plot.title = element_text(hjust = 0.5)) + geom_hline(yintercept=0,linetype="dashed",color="red") + scale_color_identity() + xlab("Birth Cohort") + ylab("Deviation from Grand Mean") + ylim(-31,41)

# Figure A65
appendixageoutpartycombinedie.rep <- data.frame(age=c(20,25,30,35,40,45,50,55,60,65,70,75,80,85,90),
                                                effect=c(1.921,0.654,0.943,0.384,-0.108,0.544,0.700,-0.283,0.63,-0.736,-1.367,-1.165,-0.474,0.743,-2.387),
                                                se=c(0.938,0.814,0.724,0.654,0.614,0.594,0.576,0.592,0.633,0.703,0.811,0.942,1.152,1.635,2.574))
appendixageoutpartycombinedie.rep$lci <- appendixageoutpartycombinedie.rep$effect - (1.96*appendixageoutpartycombinedie.rep$se)
appendixageoutpartycombinedie.rep$uci <- appendixageoutpartycombinedie.rep$effect + (1.96*appendixageoutpartycombinedie.rep$se)
appendixageoutpartycombinedie.rep$color <- ifelse(appendixageoutpartycombinedie.rep$lci<0 & appendixageoutpartycombinedie.rep$uci>0,"Gray50","Black")
appendixageoutpartycombinedie.rep$group <- 1

appendixageoutpartycombinedapci.rep <- data.frame(age=c(25,30,35,40,45,50,55,60,65,70,75,80,85),
                                                  effect=c(summary(outpartycombined.rep.raw$model)$coefficients[2:14,1]),
                                                  se=c(summary(outpartycombined.rep.raw$model)$coefficients[2:14,2]))
appendixageoutpartycombinedapci.rep$lci <- appendixageoutpartycombinedapci.rep$effect - (1.96*appendixageoutpartycombinedapci.rep$se)
appendixageoutpartycombinedapci.rep$uci <- appendixageoutpartycombinedapci.rep$effect + (1.96*appendixageoutpartycombinedapci.rep$se)
appendixageoutpartycombinedapci.rep$color <- ifelse(appendixageoutpartycombinedapci.rep$lci<0 & appendixageoutpartycombinedapci.rep$uci>0,"Gray50","Black")
appendixageoutpartycombinedapci.rep$group <- 1

figurea65.a <- ggplot(appendixageoutpartycombinedapci.rep,aes(x=age,y=effect,color=color,group=group)) + geom_point() + geom_errorbar(aes(ymin=lci,ymax=uci),size=0.5,width=0) + theme_classic() + theme(plot.title = element_text(hjust = 0.5)) + ggtitle("APCI") + geom_hline(yintercept=0,linetype="dashed",color="red") + scale_color_identity() + xlab("Age") + ylab("Deviation from Grand Mean") + ylim(-31,41)
figurea65.b <- ggplot(appendixageoutpartycombinedie.rep,aes(x=age,y=effect,color=color,group=group)) + geom_point() + geom_errorbar(aes(ymin=lci,ymax=uci),size=0.5,width=0) + theme_classic() + theme(plot.title = element_text(hjust = 0.5)) + ggtitle("IE") + geom_hline(yintercept=0,linetype="dashed",color="red") + scale_color_identity() + xlab("Age") + ylab("Deviation from Grand Mean") + ylim(-31,41)
ggarrange(figurea65.a,figurea65.b,nrow=1)

# Figure A66
appendixperiodoutpartycombinedie.rep <- data.frame(period=c(1960,1965,1970,1975,1980,1985,1990,1995,2000,2005,2010,2015),
                                                   effect=c(7.943,7.912,12.97,5.355,3.76,4.122,-2.098,-1.548,-1.83,-3.65,-14.236,-18.7),
                                                   se=c(1,0.734,0.689,0.563,0.577,0.442,0.492,0.499,0.893,0.837,0.682,0.803))
appendixperiodoutpartycombinedie.rep$lci <- appendixperiodoutpartycombinedie.rep$effect - (1.96*appendixperiodoutpartycombinedie.rep$se)
appendixperiodoutpartycombinedie.rep$uci <- appendixperiodoutpartycombinedie.rep$effect + (1.96*appendixperiodoutpartycombinedie.rep$se)
appendixperiodoutpartycombinedie.rep$color <- ifelse(appendixperiodoutpartycombinedie.rep$lci<0 & appendixperiodoutpartycombinedie.rep$uci>0,"Gray50","Black")
appendixperiodoutpartycombinedie.rep$group <- 1

appendixperiodoutpartycombinedapci.rep <- data.frame(period=c(1965,1970,1975,1980,1985,1990,1995,2000,2005,2010,2015),
                                                     effect=c(summary(outpartycombined.rep.raw$model)$coefficients[15:25,1]),
                                                     se=c(summary(outpartycombined.rep.raw$model)$coefficients[15:25,2]))
appendixperiodoutpartycombinedapci.rep$lci <- appendixperiodoutpartycombinedapci.rep$effect - (1.96*appendixperiodoutpartycombinedapci.rep$se)
appendixperiodoutpartycombinedapci.rep$uci <- appendixperiodoutpartycombinedapci.rep$effect + (1.96*appendixperiodoutpartycombinedapci.rep$se)
appendixperiodoutpartycombinedapci.rep$color <- ifelse(appendixperiodoutpartycombinedapci.rep$lci<0 & appendixperiodoutpartycombinedapci.rep$uci>0,"Gray50","Black")
appendixperiodoutpartycombinedapci.rep$group <- 1

figurea66.a <- ggplot(appendixperiodoutpartycombinedapci.rep,aes(x=period,y=effect,color=color,group=group)) + geom_point() + geom_errorbar(aes(ymin=lci,ymax=uci),size=0.5,width=0) + theme_classic() + theme(plot.title = element_text(hjust = 0.5)) + ggtitle("APCI") + geom_hline(yintercept=0,linetype="dashed",color="red") + scale_color_identity() + xlab("Year") + ylab("Deviation from Grand Mean") + ylim(-31,41)
figurea66.b <- ggplot(appendixperiodoutpartycombinedie.rep,aes(x=period,y=effect,color=color,group=group)) + geom_point() + geom_errorbar(aes(ymin=lci,ymax=uci),size=0.5,width=0) + theme_classic() + theme(plot.title = element_text(hjust = 0.5)) + ggtitle("IE") + geom_hline(yintercept=0,linetype="dashed",color="red") + scale_color_identity() + xlab("Year") + ylab("Deviation from Grand Mean") + ylim(-31,41)
ggarrange(figurea66.a,figurea66.b,nrow=1)

# Figure A67
appendixcohortoutpartycombinedie.rep <- data.frame(cohort=c(1870,1875,1880,1885,1890,1895,1900,1905,1910,1915,1920,1925,1930,1935,1940,1945,1950,1955,1960,1965,1970,1975,1980,1985,1990,1995),
                                                   effect=c(-1.343,6.423,3.699,0.185,5.446,-1.275,1.979,-0.487,1.599,1.574,0.263,-0.962,-3.063,-2.367,-3.279,-3.544,-3.63,-5.101,-2.564,-2.340,0.02,-0.963,0.213,2.156,1.768,5.593),
                                                   se=c(10.289,8.068,4.785,2.861,2.274,1.942,1.782,1.592,1.427,1.299,1.19,1.054,0.977,0.86,0.75,0.656,0.597,0.547,0.559,0.649,0.792,0.965,1.159,1.194,1.38,2.123))
appendixcohortoutpartycombinedie.rep$lci <- appendixcohortoutpartycombinedie.rep$effect - (1.96*appendixcohortoutpartycombinedie.rep$se)
appendixcohortoutpartycombinedie.rep$uci <- appendixcohortoutpartycombinedie.rep$effect + (1.96*appendixcohortoutpartycombinedie.rep$se)
appendixcohortoutpartycombinedie.rep$color <- ifelse(appendixcohortoutpartycombinedie.rep$lci<0 & appendixcohortoutpartycombinedie.rep$uci>0,"Gray50","Black")
appendixcohortoutpartycombinedie.rep$group <- 1

appendixcohortoutpartycombinedapci.rep <- data.frame(cohort=c(1875,1880,1885,1890,1895,1900,1905,1910,1915,1920,1925,1930,1935,1940,1945,1950,1955,1960,1965,1970,1975,1980,1985,1990,1995),
                                                     effect=c(as.numeric(outpartycombined.rep.raw.cohort$cohort_average$cohort_average)),
                                                     se=c(as.numeric(outpartycombined.rep.raw.cohort$cohort_average$cohort_average_se)))
appendixcohortoutpartycombinedapci.rep$lci <- appendixcohortoutpartycombinedapci.rep$effect - (1.96*appendixcohortoutpartycombinedapci.rep$se)
appendixcohortoutpartycombinedapci.rep$uci <- appendixcohortoutpartycombinedapci.rep$effect + (1.96*appendixcohortoutpartycombinedapci.rep$se)
appendixcohortoutpartycombinedapci.rep$color <- ifelse(appendixcohortoutpartycombinedapci.rep$lci<0 & appendixcohortoutpartycombinedapci.rep$uci>0,"Gray50","Black")
appendixcohortoutpartycombinedapci.rep$group <- 1

figurea67.a <- ggplot(appendixcohortoutpartycombinedapci.rep,aes(x=cohort,y=effect,color=color,group=group)) + geom_point() + geom_errorbar(aes(ymin=lci,ymax=uci),size=0.5,width=0) + theme_classic() + theme(plot.title = element_text(hjust = 0.5)) + ggtitle("APCI") + geom_hline(yintercept=0,linetype="dashed",color="red") + scale_color_identity() + xlab("Birth Cohort") + ylab("Deviation from Grand Mean") + ylim(-31,41)
figurea67.b <- ggplot(appendixcohortoutpartycombinedie.rep,aes(x=cohort,y=effect,color=color,group=group)) + geom_point() + geom_errorbar(aes(ymin=lci,ymax=uci),size=0.5,width=0) + theme_classic() + theme(plot.title = element_text(hjust = 0.5)) + ggtitle("IE") + geom_hline(yintercept=0,linetype="dashed",color="red") + scale_color_identity() + xlab("Birth Cohort") + ylab("Deviation from Grand Mean") + ylim(-31,41)
ggarrange(figurea67.a,figurea67.b,nrow=1)

# Figure A68
appendixcohortoutpartycombinedapcislopes.rep <- data.frame(cohort=c(1875,1880,1885,1890,1895,1900,1905,1910,1915,1920,1925,1930,1935,1940,1945,1950,1955,1960,1965,1970,1975,1980,1985,1990,1995),
                                                           effect=c(as.numeric(outpartycombined.rep.raw.cohort$cohort_slope$cohort_slope)),
                                                           se=c(as.numeric(outpartycombined.rep.raw.cohort$cohort_slope$cohort_slope_se)))
appendixcohortoutpartycombinedapcislopes.rep$lci <- appendixcohortoutpartycombinedapcislopes.rep$effect - (1.96*appendixcohortoutpartycombinedapcislopes.rep$se)
appendixcohortoutpartycombinedapcislopes.rep$uci <- appendixcohortoutpartycombinedapcislopes.rep$effect + (1.96*appendixcohortoutpartycombinedapcislopes.rep$se)
appendixcohortoutpartycombinedapcislopes.rep$color <- ifelse(appendixcohortoutpartycombinedapcislopes.rep$lci<0 & appendixcohortoutpartycombinedapcislopes.rep$uci>0,"Gray50","Black")
appendixcohortoutpartycombinedapcislopes.rep$group <- 1

ggplot(appendixcohortoutpartycombinedapcislopes.rep,aes(x=cohort,y=effect,color=color,group=group)) + geom_point() + geom_errorbar(aes(ymin=lci,ymax=uci),size=0.5,width=0) + theme_classic() + theme(plot.title = element_text(hjust = 0.5)) + geom_hline(yintercept=0,linetype="dashed",color="red") + scale_color_identity() + xlab("Birth Cohort") + ylab("Deviation from Grand Mean") + ylim(-31,41)

# Table A39
## Model 1
ideosort.all.raw <- APCI::temp_model(data=anes.cumulative.full,outcome='ideosort',age='agegroup',period='yeargroup',cohort='cohort',covariate=c(),family='gaussian')
## Model 2
ideosort.all.controls <- APCI::temp_model(data=anes.cumulative.full,outcome='ideosort',age='agegroup',period='yeargroup',cohort='cohort',covariate=c('gender','black','education','ideo','attendance','south'),family='gaussian')
## Model 3
ideosort.dem.raw <- APCI::temp_model(data=anes.cumulative.dem,outcome='ideosort',age='agegroup',period='yeargroup',cohort='cohort',covariate=c(),family='gaussian')
## Model 4
ideosort.dem.controls <- APCI::temp_model(data=anes.cumulative.dem,outcome='ideosort',age='agegroup',period='yeargroup',cohort='cohort',covariate=c('gender','black','education','ideo','attendance','south'),family='gaussian')
## Model 5
ideosort.rep.raw <- APCI::temp_model(data=subset(anes.cumulative.rep,agegroup!=90),outcome='ideosort',age='agegroup',period='yeargroup',cohort='cohort',covariate=c(),family='gaussian')
## Model 6
ideosort.rep.controls <- APCI::temp_model(data=subset(anes.cumulative.rep,agegroup!=90),outcome='ideosort',age='agegroup',period='yeargroup',cohort='cohort',covariate=c('gender','black','education','ideo','attendance','south'),family='gaussian')

# Figure A69
appendixideosortingieleft <- data.frame(age=c(20,25,30,35,40,45,50,55,60,65,70,75,80,85,90),
                                        effect=c(0.008,0.020,0.015,0.014,0.010,-0.003,-0.002,-0.002,-0.010,-0.016,-0.016,-0.019,-0.003,0.003,0),
                                        se=c(0.007,0.006,0.005,0.005,0.005,0.005,0.005,0.005,0.005,0.006,0.007,0.008,0.010,0.014,0.024))
appendixideosortingieleft$lci <- appendixideosortingieleft$effect - (1.96*appendixideosortingieleft$se)
appendixideosortingieleft$uci <- appendixideosortingieleft$effect + (1.96*appendixideosortingieleft$se)
appendixideosortingieleft$color <- ifelse(appendixideosortingieleft$lci<0 & appendixideosortingieleft$uci>0,"Gray50","Black")
appendixideosortingieleft$group <- 1

appendixideosortingieright <- data.frame(period=c(1975,1980,1985,1990,1995,2000,2005,2010,2015),
                                         effect=c(-0.019,-0.022,-0.029,-0.008,0.003,0.022,-0.002,0.011,0.044),
                                         se=c(0.005,0.005,0.004,0.004,0.004,0.004,0.005,0.004,0.005))
appendixideosortingieright$lci <- appendixideosortingieright$effect - (1.96*appendixideosortingieright$se)
appendixideosortingieright$uci <- appendixideosortingieright$effect + (1.96*appendixideosortingieright$se)
appendixideosortingieright$color <- ifelse(appendixideosortingieright$lci<0 & appendixideosortingieright$uci>0,"Gray50","Black")
appendixideosortingieright$group <- 1

figurea69.a <- ggplot(appendixideosortingieleft,aes(x=age,y=effect,color=color,group=group)) + geom_point() + geom_errorbar(aes(ymin=lci,ymax=uci),size=0.5,width=0) + theme_classic() + theme(plot.title = element_text(hjust = 0.5)) + ggtitle("Age") + geom_hline(yintercept=0,linetype="dashed",color="red") + scale_color_identity() + ylab("Deviation from Grand Mean") + xlab("Age") + ylim(-0.12,0.07)
figurea69.b <- ggplot(appendixideosortingieright,aes(x=period,y=effect,color=color,group=group)) + geom_point() + geom_errorbar(aes(ymin=lci,ymax=uci),size=0.5,width=0) + theme_classic() + theme(plot.title = element_text(hjust = 0.5)) + ggtitle("Period") + geom_hline(yintercept=0,linetype="dashed",color="red") + scale_color_identity() + ylab("Deviation from Grand Mean") + xlab("Year") + ylim(-0.12,0.07)
ggarrange(figurea69.a,figurea69.b,nrow=1)

# Figure A70
appendixideosortingieleft.dem <- data.frame(age=c(20,25,30,35,40,45,50,55,60,65,70,75,80,85,90),
                                            effect=c(0.043,0.040,0.020,0.013,0.001,-0.003,-0.011,-0.017,-0.029,-0.024,-0.024,-0.006,-0.024,0.010,-0.041),
                                            se=c(0.010,0.009,0.007,0.007,0.007,0.007,0.007,0.008,0.009,0.010,0.012,0.015,0.007,0.021,0.037))
appendixideosortingieleft.dem$lci <- appendixideosortingieleft.dem$effect - (1.96*appendixideosortingieleft.dem$se)
appendixideosortingieleft.dem$uci <- appendixideosortingieleft.dem$effect + (1.96*appendixideosortingieleft.dem$se)
appendixideosortingieleft.dem$color <- ifelse(appendixideosortingieleft.dem$lci<0 & appendixideosortingieleft.dem$uci>0,"Gray50","Black")
appendixideosortingieleft.dem$group <- 1

appendixideosortingieright.dem <- data.frame(period=c(1975,1980,1985,1990,1995,2000,2005,2010,2015),
                                             effect=c(-0.015,-0.022,-0.032,-0.009,-0.001,0.022,-0.006,0.010,0.053),
                                             se=c(0.007,0.007,0.005,0.005,0.005,0.006,0.007,0.006,0.008))
appendixideosortingieright.dem$lci <- appendixideosortingieright.dem$effect - (1.96*appendixideosortingieright.dem$se)
appendixideosortingieright.dem$uci <- appendixideosortingieright.dem$effect + (1.96*appendixideosortingieright.dem$se)
appendixideosortingieright.dem$color <- ifelse(appendixideosortingieright.dem$lci<0 & appendixideosortingieright.dem$uci>0,"Gray50","Black")
appendixideosortingieright.dem$group <- 1

figurea70.a <- ggplot(appendixideosortingieleft.dem,aes(x=age,y=effect,color=color,group=group)) + geom_point() + geom_errorbar(aes(ymin=lci,ymax=uci),size=0.5,width=0) + theme_classic() + theme(plot.title = element_text(hjust = 0.5)) + ggtitle("Age") + geom_hline(yintercept=0,linetype="dashed",color="red") + scale_color_identity() + ylab("Deviation from Grand Mean") + xlab("Age") + ylim(-0.12,0.07)
figurea70.b <- ggplot(appendixideosortingieright.dem,aes(x=period,y=effect,color=color,group=group)) + geom_point() + geom_errorbar(aes(ymin=lci,ymax=uci),size=0.5,width=0) + theme_classic() + theme(plot.title = element_text(hjust = 0.5)) + ggtitle("Period") + geom_hline(yintercept=0,linetype="dashed",color="red") + scale_color_identity() + ylab("Deviation from Grand Mean") + xlab("Year") + ylim(-0.12,0.07)
ggarrange(figurea70.a,figurea70.b,nrow=1)

# Figure A71
appendixideosortingieleft.rep <- data.frame(age=c(20,25,30,35,40,45,50,55,60,65,70,75,80,85,90),
                                            effect=c(-0.030,0.004,0.005,0.010,0.011,0.002,0.009,0.012,0.002,-0.005,-0.005,-0.017,-0.013,-0.011,0.020),
                                            se=c(0.009,0.008,0.006,0.006,0.006,0.006,0.006,0.006,0.006,0.008,0.008,0.009,0.011,0.016,0.027))
appendixideosortingieleft.rep$lci <- appendixideosortingieleft.rep$effect - (1.96*appendixideosortingieleft.rep$se)
appendixideosortingieleft.rep$uci <- appendixideosortingieleft.rep$effect + (1.96*appendixideosortingieleft.rep$se)
appendixideosortingieleft.rep$color <- ifelse(appendixideosortingieleft.rep$lci<0 & appendixideosortingieleft.rep$uci>0,"Gray50","Black")
appendixideosortingieleft.rep$group <- 1

appendixideosortingieright.rep <- data.frame(period=c(1975,1980,1985,1990,1995,2000,2005,2010,2015),
                                             effect=c(-0.016,-0.023,-0.027,-0.012,0.007,0.009,0.020,0.020,0.023),
                                             se=c(0.006,0.006,0.004,0.004,0.004,0.005,0.007,0.005,0.006))
appendixideosortingieright.rep$lci <- appendixideosortingieright.rep$effect - (1.96*appendixideosortingieright.rep$se)
appendixideosortingieright.rep$uci <- appendixideosortingieright.rep$effect + (1.96*appendixideosortingieright.rep$se)
appendixideosortingieright.rep$color <- ifelse(appendixideosortingieright.rep$lci<0 & appendixideosortingieright.rep$uci>0,"Gray50","Black")
appendixideosortingieright.rep$group <- 1

figurea71.a <- ggplot(appendixideosortingieleft.rep,aes(x=age,y=effect,color=color,group=group)) + geom_point() + geom_errorbar(aes(ymin=lci,ymax=uci),size=0.5,width=0) + theme_classic() + theme(plot.title = element_text(hjust = 0.5)) + ggtitle("Age") + geom_hline(yintercept=0,linetype="dashed",color="red") + scale_color_identity() + ylab("Deviation from Grand Mean") + xlab("Age") + ylim(-0.12,0.073)
figurea71.b <- ggplot(appendixideosortingieright.rep,aes(x=period,y=effect,color=color,group=group)) + geom_point() + geom_errorbar(aes(ymin=lci,ymax=uci),size=0.5,width=0) + theme_classic() + theme(plot.title = element_text(hjust = 0.5)) + ggtitle("Period") + geom_hline(yintercept=0,linetype="dashed",color="red") + scale_color_identity() + ylab("Deviation from Grand Mean") + xlab("Year") + ylim(-0.12,0.073)
ggarrange(figurea71.a,figurea71.b,nrow=1)

# Table A41
## Model 1
socialsort.all.raw <- APCI::temp_model(data=anes.cumulative.full,outcome='socialsort',age='agegroup',period='yeargroup.old',cohort='cohort.combined',covariate=c(),family='gaussian')
## Model 2
socialsort.all.controls <- APCI::temp_model(data=anes.cumulative.full,outcome='socialsort',age='agegroup',period='yeargroup.old',cohort='cohort.combined',covariate=c('gender','black','education','ideo','attendance','south'),family='gaussian')
## Model 3
socialsort.dem.raw <- APCI::temp_model(data=subset(anes.cumulative.dem,agegroup!="90"),outcome='socialsort',age='agegroup',period='yeargroup.old',cohort='cohort.combined',covariate=c(),family='gaussian')
## Model 4
socialsort.dem.controls <- APCI::temp_model(data=subset(anes.cumulative.dem,agegroup!="90"),outcome='socialsort',age='agegroup',period='yeargroup.old',cohort='cohort.combined',covariate=c('gender','black','education','ideo','attendance','south'),family='gaussian')
## Model 5
socialsort.rep.raw <- APCI::temp_model(data=subset(anes.cumulative.rep,agegroup!="90"),outcome='socialsort',age='agegroup',period='yeargroup.old',cohort='cohort.combined',covariate=c(),family='gaussian')
## Model 6
socialsort.rep.controls <- APCI::temp_model(data=subset(anes.cumulative.rep,agegroup!="90"),outcome='socialsort',age='agegroup',period='yeargroup.old',cohort='cohort.combined',covariate=c('gender','black','education','ideo','attendance','south'),family='gaussian')

# Figure A72
appendixsocialsortingieleft <- data.frame(age=c(20,25,30,35,40,45,50,55,60,65,70,75,80,85,90),
                                          effect=c(-0.021,-0.011,-0.018,-0.008,0,0,-0.008,0.003,0.002,0.005,0.012,0,0.010,0.011,0.024),
                                          se=c(0.007,0.006,0.005,0.004,0.004,0.004,0.004,0.004,0.004,0.005,0.006,0.007,0.011,0.011,0.019))
appendixsocialsortingieleft$lci <- appendixsocialsortingieleft$effect - (1.96*appendixsocialsortingieleft$se)
appendixsocialsortingieleft$uci <- appendixsocialsortingieleft$effect + (1.96*appendixsocialsortingieleft$se)
appendixsocialsortingieleft$color <- ifelse(appendixsocialsortingieleft$lci<0 & appendixsocialsortingieleft$uci>0,"Gray50","Black")
appendixsocialsortingieleft$group <- 1

appendixsocialsortingieright <- data.frame(period=c(1970,1975,1980,1985,1990,1995,2000,2005,2010,2015),
                                           effect=c(-0.022,-0.019,-0.012,-0.014,0,-0.002,0.003,0.022,0.026,0.018),
                                           se=c(0.005,0.003,0.004,0.003,0.003,0.003,0.004,0.004,0.004,0.005))
appendixsocialsortingieright$lci <- appendixsocialsortingieright$effect - (1.96*appendixsocialsortingieright$se)
appendixsocialsortingieright$uci <- appendixsocialsortingieright$effect + (1.96*appendixsocialsortingieright$se)
appendixsocialsortingieright$color <- ifelse(appendixsocialsortingieright$lci<0 & appendixsocialsortingieright$uci>0,"Gray50","Black")
appendixsocialsortingieright$group <- 1

figurea72.a <- ggplot(appendixsocialsortingieleft,aes(x=age,y=effect,color=color,group=group)) + geom_point() + geom_errorbar(aes(ymin=lci,ymax=uci),size=0.5,width=0) + theme_classic() + theme(plot.title = element_text(hjust = 0.5)) + ggtitle("Age") + geom_hline(yintercept=0,linetype="dashed",color="red") + scale_color_identity() + xlab("Age") + ylab("Deviation from Grand Mean") + ylim(-0.07,0.13)
figurea72.b <- ggplot(appendixsocialsortingieright,aes(x=period,y=effect,color=color,group=group)) + geom_point() + geom_errorbar(aes(ymin=lci,ymax=uci),size=0.5,width=0) + theme_classic() + theme(plot.title = element_text(hjust = 0.5)) + ggtitle("Period") + geom_hline(yintercept=0,linetype="dashed",color="red") + scale_color_identity() + xlab("Year") + ylab("Deviation from Grand Mean") + ylim(-0.07,0.13)
ggarrange(figurea72.a,figurea72.b,nrow=1)

## Figure A73
appendixsocialsortingieleft.dem <- data.frame(age=c(20,25,30,35,40,45,50,55,60,65,70,75,80,85),
                                              effect=c(0.001,0.007,-0.004,-0.003,0,0.008,-0.001,0.009,0,-0.002,0.002,-0.021,-0.007,0.011),
                                              se=c(0.006,0.005,0.005,0.005,0.004,0.004,0.004,0.005,0.005,0.005,0.006,0.007,0.009,0.014))
appendixsocialsortingieleft.dem$lci <- appendixsocialsortingieleft.dem$effect - (1.96*appendixsocialsortingieleft.dem$se)
appendixsocialsortingieleft.dem$uci <- appendixsocialsortingieleft.dem$effect + (1.96*appendixsocialsortingieleft.dem$se)
appendixsocialsortingieleft.dem$color <- ifelse(appendixsocialsortingieleft.dem$lci<0 & appendixsocialsortingieleft.dem$uci>0,"Gray50","Black")
appendixsocialsortingieleft.dem$group <- 1

appendixsocialsortingieright.dem <- data.frame(period=c(1970,1975,1980,1985,1990,1995,2000,2005,2010,2015),
                                               effect=c(-0.038,-0.028,-0.015,-0.014,-0.001,-0.009,-0.005,0.025,0.044,0.041),
                                               se=c(0.005,0.005,0.005,0.004,0.004,0.004,0.005,0.005,0.004,0.005))
appendixsocialsortingieright.dem$lci <- appendixsocialsortingieright.dem$effect - (1.96*appendixsocialsortingieright.dem$se)
appendixsocialsortingieright.dem$uci <- appendixsocialsortingieright.dem$effect + (1.96*appendixsocialsortingieright.dem$se)
appendixsocialsortingieright.dem$color <- ifelse(appendixsocialsortingieright.dem$lci<0 & appendixsocialsortingieright.dem$uci>0,"Gray50","Black")
appendixsocialsortingieright.dem$group <- 1

figurea73.a <- ggplot(appendixsocialsortingieleft.dem,aes(x=age,y=effect,color=color,group=group)) + geom_point() + geom_errorbar(aes(ymin=lci,ymax=uci),size=0.5,width=0) + theme_classic() + theme(plot.title = element_text(hjust = 0.5)) + ggtitle("Age") + geom_hline(yintercept=0,linetype="dashed",color="red") + scale_color_identity() + xlab("Age") + ylab("Deviation from Grand Mean") + ylim(-0.07,0.13)
figurea73.b <- ggplot(appendixsocialsortingieright.dem,aes(x=period,y=effect,color=color,group=group)) + geom_point() + geom_errorbar(aes(ymin=lci,ymax=uci),size=0.5,width=0) + theme_classic() + theme(plot.title = element_text(hjust = 0.5)) + ggtitle("Period") + geom_hline(yintercept=0,linetype="dashed",color="red") + scale_color_identity() + xlab("Year") + ylab("Deviation from Grand Mean") + ylim(-0.07,0.13)
ggarrange(figurea73.a,figurea73.b,nrow=1)

## Figure A74
appendixsocialsortingieleft.rep <- data.frame(age=c(20,25,30,35,40,45,50,55,60,65,70,75,80,85,90),
                                              effect=c(-0.046,-0.033,-0.031,-0.012,0,-0.011,-0.017,-0.008,0.002,0.009,0.021,0.020,0.021,0.015,0.07),
                                              se=c(0.008,0.007,0.006,0.006,0.006,0.005,0.005,0.005,0.006,0.007,0.008,0.009,0.011,0.015,0.026))
appendixsocialsortingieleft.rep$lci <- appendixsocialsortingieleft.rep$effect - (1.96*appendixsocialsortingieleft.rep$se)
appendixsocialsortingieleft.rep$uci <- appendixsocialsortingieleft.rep$effect + (1.96*appendixsocialsortingieleft.rep$se)
appendixsocialsortingieleft.rep$color <- ifelse(appendixsocialsortingieleft.rep$lci<0 & appendixsocialsortingieleft.rep$uci>0,"Gray50","Black")
appendixsocialsortingieleft.rep$group <- 1

appendixsocialsortingieright.rep <- data.frame(period=c(1970,1975,1980,1985,1990,1995,2000,2005,2010,2015),
                                               effect=c(-0.008,-0.014,-0.009,-0.012,0.005,0.008,0.019,0.008,0,0.004),
                                               se=c(0.006,0.005,0.005,0.004,0.004,0.004,0.005,0.007,0.005,0.006))
appendixsocialsortingieright.rep$lci <- appendixsocialsortingieright.rep$effect - (1.96*appendixsocialsortingieright.rep$se)
appendixsocialsortingieright.rep$uci <- appendixsocialsortingieright.rep$effect + (1.96*appendixsocialsortingieright.rep$se)
appendixsocialsortingieright.rep$color <- ifelse(appendixsocialsortingieright.rep$lci<0 & appendixsocialsortingieright.rep$uci>0,"Gray50","Black")
appendixsocialsortingieright.rep$group <- 1

figurea74.a <- ggplot(appendixsocialsortingieleft.rep,aes(x=age,y=effect,color=color,group=group)) + geom_point() + geom_errorbar(aes(ymin=lci,ymax=uci),size=0.5,width=0) + theme_classic() + theme(plot.title = element_text(hjust = 0.5)) + ggtitle("Age") + geom_hline(yintercept=0,linetype="dashed",color="red") + scale_color_identity() + xlab("Age") + ylab("Deviation from Grand Mean") + ylim(-0.07,0.13)
figurea74.b <- ggplot(appendixsocialsortingieright.rep,aes(x=period,y=effect,color=color,group=group)) + geom_point() + geom_errorbar(aes(ymin=lci,ymax=uci),size=0.5,width=0) + theme_classic() + theme(plot.title = element_text(hjust = 0.5)) + ggtitle("Period") + geom_hline(yintercept=0,linetype="dashed",color="red") + scale_color_identity() + xlab("Year") + ylab("Deviation from Grand Mean") + ylim(-0.07,0.13)
ggarrange(figurea74.a,figurea74.b,nrow=1)

# Table A43
## Model 1
pidstrength.all.raw <- APCI::temp_model(data=anes.cumulative.full,outcome='strong',age='agegroup',period='yeargroup.likes',cohort='cohort.likes.full',covariate=c(),family='binomial')
## Model 2
pidstrength.all.controls <- APCI::temp_model(data=anes.cumulative.full,outcome='strong',age='agegroup',period='yeargroup.likes',cohort='cohort.likes.full',covariate=c('gender','black','education','ideo','attendance','south'),family='binomial')
## Model 3
pidstrength.dem.raw <- APCI::temp_model(data=subset(anes.cumulative.dem,agegroup!="90"),outcome='strong',age='agegroup',period='yeargroup.likes',cohort='cohort.likes.full',covariate=c(),family='binomial')
## Model 4
pidstrength.dem.controls <- APCI::temp_model(data=subset(anes.cumulative.dem,agegroup!="90"),outcome='strong',age='agegroup',period='yeargroup.likes',cohort='cohort.likes.full',covariate=c('gender','black','education','ideo','attendance','south'),family='binomial')
## Model 5
pidstrength.rep.raw <- APCI::temp_model(data=anes.cumulative.rep,outcome='strong',age='agegroup',period='yeargroup.likes',cohort='cohort.likes.full',covariate=c(),family='binomial')
## Model 6
pidstrength.rep.controls <- APCI::temp_model(data=subset(anes.cumulative.rep,agegroup!="90"),outcome='strong',age='agegroup',period='yeargroup.likes',cohort='cohort.likes.full',covariate=c('gender','black','education','ideo','attendance','south'),family='binomial')

# Figure A75
appendixpidstrengthieleft.all <- data.frame(age=c(20,25,30,35,40,45,50,55,60,65,70,75,80,85,90),
                                            effect=c(-0.184,-0.138,-0.124,-0.092,-0.068,-0.027,-0.012,0.023,0.031,0.053,0.096,0.089,0.1,0.104,0.147),
                                            se=c(0.023,0.020,0.017,0.014,0.012,0.01,0.008,0.008,0.009,0.011,0.014,0.017,0.021,0.029,0.043))
appendixpidstrengthieleft.all$lci <- appendixpidstrengthieleft.all$effect - (1.96*appendixpidstrengthieleft.all$se)
appendixpidstrengthieleft.all$uci <- appendixpidstrengthieleft.all$effect + (1.96*appendixpidstrengthieleft.all$se)
appendixpidstrengthieleft.all$color <- ifelse(appendixpidstrengthieleft.all$lci<0 & appendixpidstrengthieleft.all$uci>0,"Gray50","Black")
appendixpidstrengthieleft.all$group <- 1

appendixpidstrengthieright.all <- data.frame(period=c(1950,1955,1960,1965,1970,1975,1980,1985,1990,1995,2000,2005,2010,2015),
                                             effect=c(0.037,0.059,0.047,-0.030,-0.050,-0.065,-0.018,-0.020,-0.027,-0.021,-0.023,0.022,0.047,0.044),
                                             se=c(0.023,0.019,0.017,0.014,0.011,0.008,0.008,0.007,0.009,0.011,0.014,0.017,0.019,0.022))
appendixpidstrengthieright.all$lci <- appendixpidstrengthieright.all$effect - (1.96*appendixpidstrengthieright.all$se)
appendixpidstrengthieright.all$uci <- appendixpidstrengthieright.all$effect + (1.96*appendixpidstrengthieright.all$se)
appendixpidstrengthieright.all$color <- ifelse(appendixpidstrengthieright.all$lci<0 & appendixpidstrengthieright.all$uci>0,"Gray50","Black")
appendixpidstrengthieright.all$group <- 1

figurea75.a1 <- ggplot(appendixpidstrengthieleft.all,aes(x=age,y=effect,color=color,group=group)) + geom_point() + geom_errorbar(aes(ymin=lci,ymax=uci),size=0.5,width=0) + theme_classic() + theme(plot.title = element_text(hjust = 0.5)) + ggtitle("Age") + geom_hline(yintercept=0,linetype="dashed",color="red") + scale_color_identity() + xlab("Age") + ylab("Deviation from Grand Mean") + ylim(-0.23,0.27)
figurea75.b1 <- ggplot(appendixpidstrengthieright.all,aes(x=period,y=effect,color=color,group=group)) + geom_point() + geom_errorbar(aes(ymin=lci,ymax=uci),size=0.5,width=0) + theme_classic() + theme(plot.title = element_text(hjust = 0.5)) + ggtitle("Period") + geom_hline(yintercept=0,linetype="dashed",color="red") + scale_color_identity() + xlab("Year") + ylab("Deviation from Grand Mean") + ylim(-0.23,0.27)
ggarrange(figurea75.a1,figurea75.b1,nrow=1)

## PID Strength, Democrats
appendixpidstrengthieleft.dem <- data.frame(age=c(20,25,30,35,40,45,50,55,60,65,70,75,80,85),
                                            effect=c(-0.180,-0.135,-0.112,-0.081,-0.069,-0.021,0.009,0.044,0.044,0.064,0.119,0.102,0.091,0.126),
                                            se=c(0.015,0.013,0.012,0.011,0.010,0.010,0.010,0.010,0.011,0.012,0.014,0.017,0.021,0.031))
appendixpidstrengthieleft.dem$lci <- appendixpidstrengthieleft.dem$effect - (1.96*appendixpidstrengthieleft.dem$se)
appendixpidstrengthieleft.dem$uci <- appendixpidstrengthieleft.dem$effect + (1.96*appendixpidstrengthieleft.dem$se)
appendixpidstrengthieleft.dem$color <- ifelse(appendixpidstrengthieleft.dem$lci<0 & appendixpidstrengthieleft.dem$uci>0,"Gray50","Black")
appendixpidstrengthieleft.dem$group <- 1

appendixpidstrengthieright.dem <- data.frame(period=c(1950,1955,1960,1965,1970,1975,1980,1985,1990,1995,2000,2005,2010,2015),
                                             effect=c(0.016,0.059,0.053,-0.031,-0.063,-0.074,-0.019,-0.011,-0.04,-0.018,-0.046,0.036,0.079,0.059),
                                             se=c(0.018,0.014,0.014,0.012,0.011,0.009,0.011,0.009,0.011,0.011,0.014,0.015,0.013,0.017))
appendixpidstrengthieright.dem$lci <- appendixpidstrengthieright.dem$effect - (1.96*appendixpidstrengthieright.dem$se)
appendixpidstrengthieright.dem$uci <- appendixpidstrengthieright.dem$effect + (1.96*appendixpidstrengthieright.dem$se)
appendixpidstrengthieright.dem$color <- ifelse(appendixpidstrengthieright.dem$lci<0 & appendixpidstrengthieright.dem$uci>0,"Gray50","Black")
appendixpidstrengthieright.dem$group <- 1

figurea76.a <- ggplot(appendixpidstrengthieleft.dem,aes(x=age,y=effect,color=color,group=group)) + geom_point() + geom_errorbar(aes(ymin=lci,ymax=uci),size=0.5,width=0) + theme_classic() + theme(plot.title = element_text(hjust = 0.5)) + ggtitle("Age") + geom_hline(yintercept=0,linetype="dashed",color="red") + scale_color_identity() + xlab("Age") + ylab("Deviation from Grand Mean") + ylim(-0.23,0.27)
figurea76.b <- ggplot(appendixpidstrengthieright.dem,aes(x=period,y=effect,color=color,group=group)) + geom_point() + geom_errorbar(aes(ymin=lci,ymax=uci),size=0.5,width=0) + theme_classic() + theme(plot.title = element_text(hjust = 0.5)) + ggtitle("Period") + geom_hline(yintercept=0,linetype="dashed",color="red") + scale_color_identity() + xlab("Year") + ylab("Deviation from Grand Mean") + ylim(-0.23,0.27)
ggarrange(figurea76.a,figurea76.b,nrow=1)

## Figure A77
appendixpidstrengthieleft.rep <- data.frame(age=c(20,25,30,35,40,45,50,55,60,65,70,75,80,85,90),
                                            effect=c(-0.163,-0.122,-0.123,-0.092,-0.049,-0.021,-0.032,0.004,0.024,0.047,0.077,0.081,0.119,0.083,0.166),
                                            se=c(0.026,0.022,0.019,0.017,0.015,0.013,0.012,0.012,0.013,0.015,0.018,0.021,0.026,0.037,0.053))
appendixpidstrengthieleft.rep$lci <- appendixpidstrengthieleft.rep$effect - (1.96*appendixpidstrengthieleft.rep$se)
appendixpidstrengthieleft.rep$uci <- appendixpidstrengthieleft.rep$effect + (1.96*appendixpidstrengthieleft.rep$se)
appendixpidstrengthieleft.rep$color <- ifelse(appendixpidstrengthieleft.rep$lci<0 & appendixpidstrengthieleft.rep$uci>0,"Gray50","Black")
appendixpidstrengthieleft.rep$group <- 1

appendixpidstrengthieright.rep <- data.frame(period=c(1950,1955,1960,1965,1970,1975,1980,1985,1990,1995,2000,2005,2010,2015),
                                             effect=c(0.058,0.054,0.019,-0.040,-0.035,-0.057,-0.017,-0.030,-0.003,-0.017,0.018,-0.004,0.009,0.045),
                                             se=c(0.027,0.021,0.021,0.017,0.015,0.012,0.013,0.01,0.012,0.014,0.018,0.023,0.021,0.024))
appendixpidstrengthieright.rep$lci <- appendixpidstrengthieright.rep$effect - (1.96*appendixpidstrengthieright.rep$se)
appendixpidstrengthieright.rep$uci <- appendixpidstrengthieright.rep$effect + (1.96*appendixpidstrengthieright.rep$se)
appendixpidstrengthieright.rep$color <- ifelse(appendixpidstrengthieright.rep$lci<0 & appendixpidstrengthieright.rep$uci>0,"Gray50","Black")
appendixpidstrengthieright.rep$group <- 1

figurea77.a <- ggplot(appendixpidstrengthieleft.rep,aes(x=age,y=effect,color=color,group=group)) + geom_point() + geom_errorbar(aes(ymin=lci,ymax=uci),size=0.5,width=0) + theme_classic() + theme(plot.title = element_text(hjust = 0.5)) + ggtitle("Age") + geom_hline(yintercept=0,linetype="dashed",color="red") + scale_color_identity() + xlab("Age") + ylab("Deviation from Grand Mean") + ylim(-0.23,0.27)
figurea77.b <- ggplot(appendixpidstrengthieright.rep,aes(x=period,y=effect,color=color,group=group)) + geom_point() + geom_errorbar(aes(ymin=lci,ymax=uci),size=0.5,width=0) + theme_classic() + theme(plot.title = element_text(hjust = 0.5)) + ggtitle("Period") + geom_hline(yintercept=0,linetype="dashed",color="red") + scale_color_identity() + xlab("Year") + ylab("Deviation from Grand Mean") + ylim(-0.23,0.27)
ggarrange(figurea77.a,figurea77.b,nrow=1)